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output  when  the  channel  outputs  are  either  monotonically  increasing  or  decrea¬ 
sing  in  time. 

Section  II  below  describes  the  same  model  as  in  Reference  1  except  that  the 
input  to  the  plant  is  the  difference  between  the  external  input  (pilof input) 
and  the  output  of-the  first  controller.  An  example,  with  a  plot  of  the  steady- 
state  covariances  of  the  errors  due  to  the  time  skew  between  controllers,  is 
shown  at  the  end  of  the  Section. 

An  extension  to  the  model  in  Section  II  is  developed  in  Section  III.  In  this 
extended  model,  the  first  channel  computes  two  outputs.  The  first  output  is 
the  input  to  the  plant  and  is  exactly  the  same  as  the  output  of  the  first  con¬ 
troller  of  the  model  in  Section  II.  The  second  output,  which  is  an  estimate 
of  the  output  of  the  second  channel,  is  used  to  calculate  the  error  due  to  the 
time  skew  between  the  two  controllers.  Like  the  model  in  Sectionll,  the  sec¬ 
ond  channel  computes  only  one  signal.  In  this  model,  the  inherent  errors 
depend  on  the  difference  of  the  second  output  of  the  first  channel  and  the 
output  of  the  second  channel.  An  example,  with  a  plot  of  the  steady-state  co- 
variances  of  the  second  output  of  the  first  channel  and  the  output  of  the  se¬ 
cond  channel,  is  presented  at  the  end  of  the  Section. 

An  algorithm  to  estimate  the  time  skew  between  two  asynchronous  systems  is 
described  in  Section  IV.  The  algorithm  is  based  on  the  model  in  Section  III. 
The  comparison  between  the  new  configuration  in  Sectionlll  (with  the  algorithm 
to  estimate  the  time  skew  in  SectionlV)  and  the  old  configuration  in  Section 
II  is  shown  in  Section  VI. 

Section  V  describes  the  application  of  the  new  model  in  asynchrounous  redun¬ 
dant  digital  flight  control  systems  and  Section  VII  contains  the  conclusions 
and  summary.  General  descriptions,  flowcharts,  user  instructions  and  listings 
for  all  the  software  in  this  report  are  shown  in  Appendices. 
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SECTION  I 


INTRODUCTION 

Reference  1  provides  the  background  Information  and  results 
on  the  asynchronous  operation  and  design  of  closed-loop  digital  flight 
control  systems  that  have  dual -redundant  digital  controllers.  In  the 
model  used,  the  digital  controllers  have  the  same  sample  rate  but 
there  is  a  fixed  time  skew,  or  offset  between  their  respective  sample 
times.  Also,  this  model  requires  that  the  same  channel  is  selected  as 
the  output  at  all  times.  This  latter  assumption  is  roughly  equivalent 
to  a  channel -voting  scheme  that  selects  the  upper  median  (for  a  four- 
channel  system)  or  the  lower  median  (for  a  three-  or  four-channel 
system)  as  the  output  when  the  channel  outputs  are  either  monotonically 
increasing  or  decreasing  in  time. 

Section  II  below  describes  the  same  model  as  in  Reference  1 
except  that  the  input  to  the  plant  is  the  difference  between  the 
external  input  (pilot  input)  and  the  output  of  the  first  controller. 

An  example,  with  a  plot  of  the  steady-state  covariances  of  the  errors 
due  to  the  time  skew  between  controllers,  is  shown  at  the  end  of  the 
Section. 

An  extension  to  the  model  in  Section  II  is  developed  in  Section 
III.  In  this  extended  model,  the  first  channel  computes  two  outputs. 

The  first  output  is  the  input  to  the  plant  and  is  exactly  the  same  as 
the  output  of  the  first  controller  of  the  model  in  Section  II.  The 
second  output,  which  is  an  estimate  of  the  output  of  the  second  channel. 
Is  used  to  calculate  the  error  due  to  the  time  skew  between  the  two 
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controllers.  Like  the  model  in  Section  II,  the  second  channel  computes 
only  one  signal.  In  this  model,  the  inherent  errors  depend  on  the 
difference  of  the  second  output  of  the  first  channel  and  the  output  of 
the  second  channel.  An  example,  with  a  plot  of  the  steady-state 
covariances  of  the  second  output  of  the  first  channel  and  the  output 
of  the  second  channel,  is  presented  at  the  end  of  the  Section. 

An  algorithm  to  estimate  the  time  skew  between  two  asynchronous 
systems  is  described  in  Section  IV.  The  algorithm  is  based  on  the 
model  in  Section  III.  The  comparison  between  the  new  configuration 
in  Section  III  (with  the  algorithm  to  estimate  the  time  skew  in  Section 
IV)  and  the  old  configuration  in  Section  II  is  shown  in  Section  VI. 

Section  V  describes  the  application  of  the  new  model  in  asyn¬ 
chronous  redundant  digital  flight  control  systems  and  Section  VII 
contains  the  conclusions  and  summary.  General  descriptions,  flow 
charts,  user  instructions  and  listings  for  all  the  software  in  this 
report  are  shown  in  Appendices. 


SECTION  II 


STATE  EQUATIONS,  COVARIANCE,  AND  EXAMPLES 
FOR  BASIC  MODEL 

The  model  illustrated  in  Figure  1  is  labelled  the  basic  model; 
the  assumptions,  techniques,  and  style  of  analysis  are  the  foundation 
for  the  new  model  described  in  Section  III  of  this  report.  The  basic 
model  is  similar  to  the  model  in  Reference  1,  except  that  the  input 
to  the  plant  is  the  difference  between  the  external  input  (pilot 
input)  and  the  output  of  the  first  controller,  while  in  the  model  of 
Reference  1,  the  input  to  the  plant  is  the  output  of  the  first  controller. 

1.  SYSTEM  CONFIGURATION  AND  THE  DYNAMIC  EQUATION 

The  system  configuration  for  this  closed-loop  dynamic  system 
consists  of  a  conti nuous-time  plant  and  dual -redundant,  single-rate 
discrete-time  controllers.  The  plant  output  is  sampled  by  each  of  the 
controllers,  using  a  common  sample  period  but  having  a  fixed  time  skew 
between  them.  The  output  of  one  of  the  controllers  serves  as  the 
piecewise-constant  input  to  the  plant,  along  with  an  external  input. 

The  plant  equations  include  aircraft,  sensor,  and  actuator  dynamics, 
as  well  as  any  dynamics  associated  with  the  pilot  input  and  wind-gust 
model  input.  The  plant  equations  are  assumed  to  be  in  the  form 

xp  =  ^pxp  +  ®pup  (2-1) 

Yp  *  CpXp  (2-2) 


where 


T  :  SAMPLE  PERIOD 
x  :  SKEW  TIME 


FIGURE  I  :  BLOCK  DIAGRAM  FOR  THE  BASIC  MODEL 
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Xp  =  plant  state  vector  (np  x  1) 
up  =  plant  input  vector  (nUp  x  1) 

Wp  =  external  input  vector  (nup  x  1) 
ycl=  controller  1  output  vector  (nup  x  1) 
yp  *  plant  output  vector  (nop  x  1) 

Ap  =  plant  state  matrix  (rtp  x  np) 

Bp  =  plant  input  matrix  (np  x  nup) 

Cp  =  plant  output  matrix  (n0p  x  np) 

The  solution  to  Equation  1  is 

t 

xp(t)  =  ♦(t,t0)xp(to)  +  f  4>(t,s)BpUp(s)ds  (2-3) 

to 

where  4»(t,tQ)  is  the  state  transition  matrix  and  for  constant  Ap 
is  given  by 

*(t,t0)  =  exp  [Ap(t-t0)3  (2-4) 

The  plant  input  up(t)  is  piecewise-constant  over  a  given  sampling 
interval;  i.e., 

up(t)  =  up(tk)  tk  <  t  <  tk+i 

and  so  for  t  ■  tk,  t  *  tk+^,  and  tk+^  -  tk  *  T,  the  second  term  in 
Equation  (2-3)  can  be  written  as 

lk+l 

f  ♦(tk+i,s)BpUp(s)ds  -  *(t|(+i,tk)up(tk)  (2-5) 

tk 
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where 


^Ctk+1  »tk)  =  /  exp  [Ap(t)]Bpdt  (2-6) 

o 

Substitution  of  (2-6)  into  (2-3)  gives 

xp(tk+l )  “  ^(^k+1  *^k^xp(t|()  +  ’/'(tk+i  *t|c)Up(t|{)  (2-7) 

for  k  =  0,  1 ,  .  .  . 

The  discrete- time  equations  for  controller  #1  are 


xcl (tk+l 

)  =  ^cxcl^tk)'+  ^cucl(^k) 

(2-8) 

yci^k) 

Hcxd(tk)  +  Ecuci(tk) 

(2-9) 

for  k  =  0,  1 ,  .  . 

• 

and  for  controller 

#2 

xc2(tk+l 

+  T)  ■  FcxC2^k  +  T)  +  Gcuc2^tk  +  T) 

(2-10) 

yc2(tk  + 

t)  =  HcxC2(tk  +  t)  +  Ecuc2 ( tk  +  t) 

(2-11) 

for  k  =  0 ,  1 ,  .  .  . 
where 

xcl  *  controller  1  state  vector  (nc  x  1) 
ycl  =  controller  1  output  vector  (nup  x  1) 
xc2  =  controller  2  state  vector  (nc  x  1) 
yC2  =  controller  2  output  vector  (nup  x  1) 
F£  *  controller  state  matrix  (nc  x  nc) 


6 


Gc  =  controller  control  input  matrix  (nc  x  nop) 

Hc  *  controller  output  matrix  (states)  (nup  x  nc) 

Ec  =  controller  output  matrix  (inputs)  (nup  x  nop) 

The  plant  (aircraft,  actuator,  and  sensor  dynamics)  and  the 


controllers  are  related  by  the  equations 

up(tk)  “  wp(tk)  -  yci(tk)  (2-12) 

Uciitic)  s  yp(tk)  (2_13) 

ucZ(tk  +  T)  =  *  T )  (2-14) 

Substitution  of  Equation  (2-12)  into  (2-7)  gives 

xp^k+l^  *  ♦(*k+l  »tk)xp(tk)  + 

4>(tk+i.tk)[wp(tk)  -  yci(tk)3  (2"15) 

and 

yp^k+l)  =  Cpxp(tk+i )  (2-16) 


The  quantity  xp(tk  +  t)  can  be  written  using  the  solution  to 
Equation  (2-7)  as 

xp(tk  +  t)  »  *(tk  +  T,tk)xp(tk)  + 

*(tk  +  T,tk)[wp(tk)  -  yci(tk)]  (2-17) 

and 

yp(tk  +  t)  -  Cpxp(tk  +  t)  (2-18) 
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The  piecewise-constant  inherent  error  e(t)  is  written  in  two 
parts,  e^(t)  and  eB(t)  as 

eA^)  =  ycl  ‘  yc2^k  +  *)  (2-19) 

for  tk  +  t  <_  t  <  t^ ,  0  <  t  <  T,  k  *  0,  1,  .  .  .  ,  and 

eB(t)  3  Yc^^k+l^  “  yC2^k  +  (2-20) 

for  tk+i  <  t  <  tk+1  +  t,  0  <  t  <  T,  k  =  0,  1,  .  .  . 

2.  COVARIANCE  ANALYSIS 

Let  the  input  wp(tk)  be  a  Gaussian  white  noise  random  process 
with  zero  mean,  which  is  independent  of  x(0)  (Reference  3).  Then, 
let  PfiA  and  PgB  be  the  covariance  of  the  errors  eA  and  eg,  respectively. 
Thus, 

PeA(t)  =  E[eA(t)  e aT ( t ) ]  (2-21) 

for  tk  +  x  <_  t  <  tk+.j ,  0  <  x  <  T,  k  =  0,  1 ,  .  .  .  ,  and 

PeB(t)  -  E[eB(t)  eBT(t)]  (2-22) 

for  *k+l  tk+l  +  t»0<t<T,  k  =  0,  1,.  .  . 

Since  the  input  wp  is  a  Gaussian  white  noise  and  the  controllers 
are  discrete  time  controllers,  the  inherent  errors  e^  and  eB  will  be 
random  variables  (g^| ,  ^A2*  •  •  •  ,  eA|^},  {egj ,  ^B2*  •  *  *  *  eg^l, 
as 

eA1  Is  the  error  In  the  interval  tk  +  t  and  tk+1 
Is  the  error  in  the  Interval  tk+i  +  t  and  tk+2 
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eAN  is  the  error  in  the  interval  t^.^  +  r  and  t^+n 
and 

egl  is  the  error  in  the  interval  t^+j  and  tk+1  +  t 

eg2  is  the  error  in  the  interval  t^+2  and  t^+2  +  t 

eBN  is  the  error  in  the  interval  t^+fj  and  t^+m  +  t 

Let  e^  and  eg  represent  the  sample  means  based  on  N  samples  of 

{6)\1 1  ®A2*  •  •  •  »  ®AN^  (fifii »  ®B2*  *  •  *  *  eBN^‘  Thus j 

FA  3  I  teAl  +  eA2  +  '  •  *  +  eAN^  (2-23) 

N 

and 

=  I  [®gi  ®g2  ■*■•••"*■  eBN]  (2-24) 

N 

From  these  two  equations  (2-23)  and  (2-24),  the  sample  covariances 
of  the  errors  based  on  the  interval  from  t^  to  t^  are 

peA  *  I  j,  C'AI  -  V  C'AI  -  *yT  <2-25> 

N  1=1 

and 

%B  "  -  ‘  t‘A)  '  *A^  t'Af  •  *fr|T 

N  1*1 

For  the  steady  state  sample  covariance  of  errors,  k  is  the  value 
when  the  system  is  in  steady-state. 


9 


3.  EXAMPLE 

As  the  example,  consider  the  closed-loop  system  shown  in  Figure 
2  with  a  second-order  plant  and  a  first-order  controller.  Assume  wp 
to  be  Gaussian  white  noise  with  zero  mean  and  variance  s  1  and  let 
the  sample  period  T  equal  0.0125  seconds. 

From  the  block  diagram  in  Figure  2,  the  plant  is  described  by 


A 


P 


s 


0 

0 


1 


-10 


V 


0 

200 


and 


The  description  of  the  digital  controllers  is  obtained  by  starting 
with  the  Laplace  transform  transfer  function  of  an  analog  controller 
and  performing  the  Tustin  transformation  (also  called  the  bilinear 
transformation)  to  obtain  the  z  transformation  and  the  discrete-time 
state  equations. 

The  continuous  controller  transfer  function  is  1  *  0«Q3s  . 

1  +  0.02s 

The  substitution  s  =  £  z~^  performs  the  Tustin  transformation. 

0.03  x  2  z-1 
1  +  T  zTT 

0.02  x  2  z-1 
T  z+1 


T  Z+1 

This  yields 

1  +  0.03s  _ 

1  +  0.02s 
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FIGURE  2 


Then,  the  transfer  function  can  be  written  as 

Q.Q4T _ 

Xcl  =  T  +  0.06  -  Tz  +  0.08T  +  0.0016 

Yp  T  +  0.04  z  .  0.04  -  T 

0.04  +  T 

where  Xc-| (z)  is  the  controller  output  and  Yp(z)  is  the  controller 

input,  which  is  also  the  plant  output.  A  block  diagram  for  digital 

controller  1  appears  in  Figure  3. 

The  state  equations  corresponding  to  Figure  3  are 


Xcl^k+l)  =  Fcxcl  +  Gcyp(tk) 
ycl(tk)  =  Hcxcl(tk)  +  Ecyp(tk) 

where  tk+]  -  tk  =  T  and 


F  *  0.04  -  T 
c  0.04  +  T 

Gc  *  0.04T 

T2  +  0.08T  +  0.0016 


Ec  ■  0.06  +  T 
0.04  +  T 

The  state  transition  matrix  ♦(tk+1,tk)  from  equation  (2-4)  is 

-10T 


♦^k+l’V 


1  TO 
0 


1  -  1 


TO 


-10T 


The  steady  state  sample  variance  of  errors  Pe^ss  and  Pegss  from 
Appendix  B  are  plotted  in  Figure  4  as  a  funtion  of  t.  The  diagrams  to 
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the  right  of  each  plot  show  the  times  at  which  the  controller  outputs 
are  sampled  for  the  calculation  of  e^  and  eg.  The  sample  variances 
are  largest  when  the  times  at  which  ycl  and  yc2  change  are  farthest 
apart,  as  expected.  The  results  indicate  that  some  combination  of 
the  two  measurements  of  the  channel  inherent  errors  may  be  less  affected 
by  the  amount  of  skews  than  either  error  taken  alone. 
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SECTION  III 


STATE  EQUATIONS,  COVARIANCE,  AND 
EXAMPLE  FOR  NEW  MODEL 

According  to  the  example  of  the  basic  model  in  Section  II,  the 
tolerance  value  for  two-channel  operation  is  the  maximum  value  of  the 
steady-state  sample  covariance  Pe»ss  and  is  reached  at  t  =  T.  If 
the  steady-state  covariance  is  greater  than  the  tolerance  value,  the 
error  is  due  to  a  failed  channel ;  if  it  is  less  than  the  tolerance, 
the  error  is  the  inherent  error  due  to  sampling  skew.  However,  the 
above  choice  for  the  tolerance  value  may  not  be  the  best  one  to 
distinguish  the  inherent  error  from  the  error  due  to  a  failed  channel. 

For  example,  if  the  time  skew  of  the  second  channel  is  small  and 
the  sample  variance  is  greater  than  the  average  variance  for  that  time 
skew,  then  it  is  likely  that  there  is  a  failed  channel.  If  that  sample 
variance,  which  is  greater  than  the  average  variance  for  that  skew 
time,  is  less  than  the  tolerance  value  (PeAss  at  t  =  T),  the  basic 
model  in  Section  II  would  indicate  no  failure. 

To  reduce  the  effect  above,  one  possibility  is  for  channel  1 
to  compute  an  approximation  to  the  current  output  of  channel  2.  The 
difference  between  the  true  output  of  channel  2  and  the  estimated 
output  of  channel  2  will  be  close  to  zero,  assuming  that  the  estimate 
is  a  good  one.  The  tolerance  value  for  this  approach  can  be  a  small 
value  and  it  is  equal  to  the  maximum  steady-state  sample  covariance  of 
the  difference  between  the  estimated  value  and  the  actual  value.  The 
model  described  In  this  section  estimates  the  output  of  channel  2 
from  the  input  of  channel  1  and  the  details  of  this  approach  are  in 
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the  following  subsections. 


1 .  SYSTEM  CONFIGURATION  AND  DYNAMIC  EQUATIONS 

The  closed-loop  system  configuration  for  this  new  model  appears 
in  Figure  5,  and  is  almost  the  same  as  the  basic  model  in  Section  II. 
The  output  yp  of  the  plant,  which  consists  of  the  aircraft,  the  sensor, 
and  the  control  actuator  dynamics,  is  sampled  by  each  of  two  digital 
controllers.  They  use  the  same  fixed  sample  period  T,  but  there  is 
a  constant  skew  t  between  the  starting  points  of  the  two  samplers. 

In  Figure  5,  the  input  of  channel  1,  ypU^),  is  used  to  compute 
yci(t|<)  and  y£2(tk  +  T*)»  an  estimate  of  yc2Uk  +  t)  (t*  is  an  estimate 
of  t).  The  difference  between  the  external  input  wp  and  yc-j  is  the 
input  to  the  plant.  The  block  named  OBSERVER  (See  Appendix  D  for 
details)  computes  XpU^),  an  estimate  of  XpU^)  based  on  the  two 
quantities  ypUj^)  and  up(t|c_]),  which  are  the  previous  inputs  of 
channel  1  and  the  plant,  respectively.  yp(tk  +  T*)  is  an  estimate  of 
yp(tk  +  t)  calculated  from  t*,  x£(tfc),  and  yc](tk);  this  estimate  is 
an  input  to  the  block  named  2nd  DIGITAL  CONTROLLER  #1.  Finally, 
y*2(tk  +  t*)  is  computed. 

As  in  the  basic  model,  the  plant  equations  include  the  aircraft, 
sensor,  and  actuator  dynamics,  as  well  as  any  dynamics  associated  with 
the  pilot  input  and  the  wind-gust  model  input.  The  plant  equations 
are  assumed  to  be  in  the  form 

Xp  =  ^pxp  +  ®pup  (3-1 ) 

yp  *  CpXp  (3-2) 
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AIRCRAFT  AND  ACTUATOR  DYNAMICS 


Without  showing  the  details  of  the  derivation  (the  details  are  avail¬ 
able  in  Section  II),  the  solution  of  equation  (3-1)  is 

xp(tj<+i)  =  <j>(tk+^,  tk)xp(tk)  +  i|»(tk+],  tk)up(tk)  (3-3) 

where 

up(tk)  =  wp(tk)  -  yci (tk)  (3-4) 

The  first  function  of  channel  1  is  to  compute  the  signal  that 
is  fed  back  to  the  plant  according  to  the  equation 

yd^k)  =  8cxcl^k)  +  ^c^p^k^  (3-5) 

and 

xc](tk+i)  =  Fcxcl^k^  +  GcyP^k)  (3-6) 

for  k  =  0,  1 ,  2,  .  .  . 

The  second  function  is  to  compute  the  signals  that  are  used  to 
calculate  the  inherent  error  according  to 

y*2(tk  +  t*)  =  Hcx*z(tk  +  t*)  +  Ecyp(tk  +  x*)  (3-7) 

and 

xc2^k+l  +  T*)  =  Fcxc2^k  +  T*)  +  6cyp(tk  +  T*)  (3-8) 

for  k  =  0 ,  1 ,  2 ,  .  .  .  Here , 

y*2(tk  +  x*)  is  the  estimate  of  yc2(tk  +  T) 

xc2^k  +  T*)  *s  the  est,,mate  oF  xc2(tk  +  T) 
and 

y*(tk  +  x*)  is  the  estimate  of  yp(tk  +  x) 
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Note:  All  *  variables  are  in  channel  1  and  the  computations  are 
done  at  the  times  tk,  tk+],  tk+2,  ....  instead  of  tk  +  t, 
lk+l  +  T»  tk+2  +  t,  .  .  . 

Channel  2  computes  the  signal  that  is  to  calculate  the  inherent 
error  according  to 

yC2^k  +  T)  =  HcxC2(fck  +  T)  +  Ecyp(*k  +  T)  (3-9) 

and 

xc2(tk+l  +  T)  =  Fcxc2^k  +  T)  +  G^p^k  +  T)  (3-10) 

From  equation  (3-3),  we  can  write  the  equation  of  xp(tk  +  t) 
as 


xp(tk  +  t)  =  *(tk  +  t,  tk)xp(tk)  +  *(tk  +  t,  tk)up(tk)  (3-11) 
and  the  input  of  controller  #2  is 


is 


From 


yp(tk  +  t)  =  Cpxp(tk  +  t)  (3-12) 

Appendix  0,  the  equation  of  the  observed  state  xp(tk+j) 


xp^k+l^  =  ^(^k+1’  ^k^xp^k^  +  ♦(^k+l*  ^k^up(^k^ 

+  *1^+1  ’  WV  <3-13> 


where 


and 


xp(tk)  is  the  estimate  of  xp(tk) 


T 

♦l (tk+1 »  tk)  =  /  EXP(Ap(t))6edt 
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9 


where  Gg  is  the  feedback  matrix  of  the  observer  which  y£  will  approach 

V 

From  equation  (3-9),  the  equation  of  x*(tk  +  t  )  can  be  written 
as 

x*(tk  +  T*)  =  ♦Uk  +  t*.  tk)x*(tk) 

+  *(tk  +  t*.  tk)up(tk)  (3-15) 

where  t*  is  the  estimate  of  t. 

In  computing  +  T*)»  the  variable  y*(tk  +  t*)  is  calculated 

from  the  equation 

yj(tk  +  t*)  -  CpXp(tk  +  t*)  (3-16) 

The  piecewise-constant  Inherent  error  e(t)  is  written  in  two 
parts,  eA(t)  and  eg(t)  as 

=  ytf(tk  +  T*)  "  yc2^k  +  T)  (3-17) 

for  tk  +  <^  t  •<  tk+^ ,  0  x  <  T,  k  =  0,  1 ,  .  .  •  ,  and 

eB^)  s  y^^k+l  +  T*)  ’  yc2^k  +  T)  (3-18) 

for-  tk+1  <  t  <  tkt,  *  0  <  t  <  T,  k  •  0,  1,  .  .  . 

2.  COVARIANCE  ANALYSIS 

As  in  the  basic  model,  Wp  is  the  Gaussian  white  noise,  which  has 
zero  mean  and  variance  *  1.  Thus,  the  samples  covariance  of  errors 
eA  and  eg  are 
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(3-19) 


%A(t>  -  J  j,  (eA1  -  •*)  (eAI  -  eA)T 
for  tk  +  t  £  t  <  tk+1 ,  0  <_  t  <  T,  k  =  0,  1 ,  .  .  .  ,  and 

PeB^)  =  1  ?  (eAi  "  eA)T  (3-20) 

N  1=1  A1  A 

for  tk+1  <  t  <  tk+]  +  x,  0  <  t  <  T,  k  =  0,  1,  .  .  . 

3.  EXAMPLE 

As  the  example,  consider  the  system  In  Figure  2  of  Section  II. 
peAss  and  peBss  e9uations  (3-19)  and  (3-20)  are  plotted  as  a  function 
of  t  for  t*  =  0,  T/5,  2T/5,  3T/5,  4T/5,  and  T  respectively  in  Figure  6 
(from  Appendix  E).  The  diagrams  to  the  right  of  each  plot  show  the 
times  at  which  the  controller  output  are  used  for  the  calculation  of 
eA  (eq.  3-17)  and  eg  (eq.  3-18).  The  sample  variance  (PeAss)  are  largest 
when  t*  is  farthest  apart  from  t  (or  y^^k  +  T*)  1S  farthest  apart 
from  yc2(tk  +  as  expected. 

From  the  plot,  if  t*  is  equal  to  t,  the  inherent  error  (eft)  of 
this  model  is  zero  and  the  main  disadvantage  of  the  asynchronous  operation 
will  be  eliminated.  If  t*  is  close  to  t,  the  inherent  error  (eA)  is 
a  small  value.  Then  the  deficiency  of  the  basic  model,  which  is  de¬ 
scribed  at  the  beginning  of  this  section,  will  be  reduced.  According  to 
the  plot,  the  sample  variance  of  eA  (peAss)  Is  directly  proportional 
to  the  difference  between  t  and  t*.  Then  t  can  be  estimated  by  comparing 
the  sample  variances  of  eA  for  the  values  of  t*  in  [0,  T]  and  t* 
which  corresponds  to  the  smallest  covariance  will  be  the  estimate  of  t. 

The  next  section  will  describe  the  detail  of  the  algorithm  for  estimating 
t  by  using  the  result  discussed  above. 
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SECTION  IV 


ALGORITHM  FOR  THE  ESTIMATOR  t* 

According  to  the  example  in  Section  III,  the  steady-state  sample 
covariance  of  eft  (eq.  3-19)  appears  to  be  directly  proportional  to  the 
difference  between  t  and  t*.  That  is,  when  the  difference  between  t 
and  t*  is  large,  the  steady-state  sample  covariance  of  e^  is  large;  when 
the  difference  is  small,  the  covariance  is  small;  and  when  the  difference 
is  zero,  the  covariance  is  zero.  This  relationship  is  the  basis  for  the 
technique,  described  in  this  section,  to  estimate  x*.  The  technique 
uses  the  model  of  Figure  7  and  the  assumption  that  the  steady-state 
sample  covariance  of  depends  on  the  difference  between  x  and  x*. 

1.  DESCRIPTION  OF  THE  ALGORITHM 

The  basic  procedure  for  estimating  x*  by  using  the  model  of  Figure 
7  is  to  change  x*  in  an  iterative  manner  until  the  smallest  covariance 
of  e^  is  obtained.  Let  the  single  variable  NT  be  the  number  of  sub¬ 
intervals  in  the  interval  (0,T)  so  that  the  length  of  each  subinterval 

in  (0,T)  is  equal  to  ^  and  (0,T)  is  divided  into  0,  T,  .  .  .  , 
nt  i 

Kl  "  1  x  T,  T.  Let  0  and  T  be  the  first  lower  and  upper  limits  in  which 
NT 

x  lies.  Let  x*  be  the  estimate  of  t.  By  comparing  the  steady-state 

sample  covariance  of  e^  when  x*  is  the  midvalue  between  the  lower  and 

the  upper  limits  and  that  of  eA  when  x*  Is  the  value  which  is  greater 

than  the  midvalue  by  ,  one  can  determine  whether  to  Increase  or  to 

NT 

decrease  x*.  to  reduce  the  steady-state  sample  covariance.  If  the 
previous  steady-state  sample  covariance  of  eA  (eq.  3-19)  when  t*  Is 
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IE  MODEL  FOR  ESTIMATING 


1 


the  midvalue  is  less  than  the  current  steady-state  sample  covariance 

of  e^,  when  x*  is  greater  than  the  midvalue  by  _J_  ,  then  t  must  be 

NT 

between  the  lower  limit  and  the  midvalue.  Therefore,  this  midvalue 

is  selected  as  the  new  upper  limit  of  the  new  interval  of  t  while  the 

lower  limit  is  unchanged.  If  the  previous  steady-state  sample  covariance 

of  e^  is  greater  than  the  current  steady-state  sample  covariance  of 

eft,  then  the  lower  limit  is  updated  with  the  midvalue  while  the  upper 

limit  is  maintained.  (Note  that  the  values  of  the  lower  limits,  the 

upper  limits,  and  the  midvalues  are  restricted  to  the  values  0,  _J_  , 

NT 

££,  .  .  .  ,  T.)  This  procedure  is  repeated  until  the  new  midvalue 
NT 

differs  from  the  previous  midvalue  by  T.  Then  the  resulting  interval 

NT 

is  the  smallest  interval  in  which  t  lies  and  any  value  of  x*  in  this 

,  smallest  interval  can  be  used  as  the  estimate  of  t. 

When  the  smallest  interval  in  which  t  lies  is  obtained,  the 

lower  limit  or  the  upper  limit  will  be  the  previous  midvalue.  Then 

the  last  current  midvalue  will  differ  from  the  lower  or  the  upper 

limits  by  However,  the  midvalue  is  not  necessarily  equal  to  the 
NT 

exact  midpoint  between  the  lower  and  the  upper  limits.  Figure  8 
shows  the  two  possible  locations  of  the  midvalue.  (Note:  the  midvalue 
in  this  report  is  the  average  value  of  the  lower  and  upper  limits  in 
which  the  average  value  is  a  truncated- Integer  division.)  Since  the 
average  value  of  the  lower  and  upper  limits  is  a  truncated- Integer 
division,  then  there  are  three  smallest  intervals  of  x  which  is  shown 
in  Figure  9. 
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FIGURE  9  :  THREE  POSSIBLE  SMALLEST  INTERVALS  OF  t 
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This  technique  is  similar  to  the  'HALF- INTERVAL  SEARCH'  which 
is  a  method  for  obtaining  an  approximate  solution  to  an  equation 
f(x)  =  0  and  it  is  available  in  almost  every  numerical  analysis  book. 

A  flowchart  of  this  procedure  appears  in  Figure  10.  Before 
discussing  the  flowchart,  let  us  define  the  variables  which  are  used 
in  this  flowchart. 

As  mentioned  before,  NT  is  the  number  of  subintervals  in  (0,T) 

and  instead  of  using  the  real  values  of  the  subintervals  in  (0,T); 

namely,  )  x  T  )  x  T  ,  .  .  .  ,  (NT1  -1 )  x  T,  it  is  more 

(NT1-1 )  (NT1-1 )  (NT1-1) 

convenient  to  refer  to  the  numbers  1,  2,  .  .  .  ,  NT1.  N4  and  N5  are 
the  integers  which  represent  the  lower  and  upper  limits  of  t  respec¬ 
tively.  The  lower  and  the  upper  limits  of  the  interval  which  t  lies. 

N3  is  the  truncated- integer  midvalue  of  N4  and  N5. 

NTAU1  and  NTAU2  are  the  integers  which  represent  t*  at  the 
updated  midvalue  and  the  previous  midvalue  respectively.  NTAU3  is 
also  the  integers  which  represents  t*  at  the  midvalue  plus  one.  PEASS1 
and  PEASS2  are  defined  to  be  the  steady-state  sample  covariances  of 
e^  which  correspond  to  NTAUl  and  NTAU3  respectively. 

The  first  step  of  the  flowchart  shows  the  initialization  of 
the  key  variables. 

The  second  step  shows  the  computation  of  the  covariance  of  eA 
when  t*  is  equal  to  the  midvalue  of  the  interval.  NTAUl  represents 
the  midvalue  and  PEASS1  is  the  corresponding  covariance. 

The  third  step  provides  the  decision  used  in  terminating  the 
routine.  This  routine  will  terminate  when  the  current  midvalue  (NTAUl) 
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NTAU2=0 , N4»0 , N5=NT1 , 
N3=(N4+N5)/2 


NTAU1-N3 ; COMPUTE  PEASS1 


INTAU1-NTAU2I  : 


N3-N3+1 , NTAU3-N3 ; 
COMPUTE  PEASS2 
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FIGURE  10  :  FLOWCHART  OF  THE  DETAILS  OF 

THE  PROCEDURE  FOR  ESTIMATING  T* 


FIGURE  10  :  (cont.) 
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is  one  subinterval  apart  from  the  previous  midvalue  (NTAU2). 

The  fourth  step  provides  another  value  of  the  covariance  when 
r*  is  equal  to  NTAU3. 

The  fifth  step  takes  care  of  the  comparison  between  the  steady- 
state  sample  covariance  obtained  from  step  2  and  that  obtained  from 
step  4.  The  decision  is  made  in  this  step  in  order  to  select  a  new 
interval  of  x.  The  steady-state  sample  covariance  of  e^  is  directly 
proportional  to  the  difference  between  t  and  x*;  therefore,  if  PEASS1 
is  greater  than  PEASS2,  then  the  lower  limit  N4  is  updated  with  the 
midvalue  while  the  upper  limit,  N5,  is  unchanged.  If  PEASS1  is  less 
than  PEASS2,  then  N5  will  be  updated  with  the  midvalue  while  N4  is 
unchanged.  The  procedure  goes  back  to  step  2  and  repeats  until  the 
condition  in  the  third  step  is  met. 

The  remainder  of  the  flowchart  (after  the  difference  of  the  current 
midvalue  and  the  previous  midvalue  is  equal  to  one)  shows  the  details 
of  the  technique  for  estimating  t.  As  discussed  at  the  beginning  of 
this  section,  the  estimate  of  x  can  be  selected  by  comparing  the  steady- 
state  sample  covariance  of  every  quantized  number  between  the  updated 
N4  and  N5.  The  integer  in  this  interval  which  corresponds  to  the 
smallest  covariance  denotes  the  estimate  of  x. 

If  NT  is  an  integer  power  of  2;  i.e., 

NT  =  2V 

then  there  is  only  one  possible  smallest  interval  in  which  x  lies 
and  it  is  type-b  smallest  interval,  which  is  shown  in  Figure  9-a. 
Furthermore,  the  maximum  number  of  iterations  to  get  the  estimate  of 
x  for  any  value  of  NT  can  be  estimated  by  the  maximun  number  of 
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*  smallest  interval  type-b (FIGURE  9) 

*  smallest  interval  type-c (FIGURE  9) 
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FIGURE  11  :  DIAGRAM  OF  THE  NUMBER  OF 
INTE RATIONS  OF  (a)  NT-8, 
(b)  NT-9,  (c)  NT- 10 
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FIGURE  11  :  (cont.)  (d)  NT-11,  (e)  NT-12 
(f)  NT-13 
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iterations  when  NT  is  equal  to  an  Integer  power  of  2.  For  example, 
if  NT  is  8,  then  the  maximum  Iterations  is  3  as  shown  in  Figure  11 -a. 

The  last  Iteration  (3rd  iteration)  is  the  smallest  interval  of  t  when  the 
difference  between  the  previous  midvalue  and  the  current  midvalue  is 
one.  All  these  intervals  are  type-b  smallest  intervals.  If  NT  is  increased 
by  1  (NT  =  9),  then  the  last  Iteration  (3rd  iteration)  contains  three 
type-b  smallest  interval  and  one  type-c  smallest  interval  and  it  is  shown 
in  Figure  11-b.  The  same  as  NT  =  8,  the  maximum  iterations  of  NT  =  9 
is  equal  to  3. 

The  situation  when  NT  is  10,  11,  and  12  are  shown  in  Figure  11 -c, 
d,  and  e,  respectively,  and  the  maximum  iterations  of  these  values  of 
NT  is  still  equal  to  3.  When  NT  is  13,  the  maximum  iterations,  which 
is  shown  in  Figure  11-f,  is  4.  From  these  examples,  the  approximate 
maximum  iterations  of  any  values  of  NT  between  the  midvalue  of  2V"^ 
and  2^  and  the  midvalue  of  2^  and  2^  is  equal  to  V. 

Another  approach  for  calculating  the  estimate  of  t  is  to  determine 
the  steady-state  sample  covariance  of  e^  for  successive  values  of 
t*  in  (0,T)  until  that  covariance  begins  to  increase.  The  t*  which 
corresponds  to  the  smallest  value  of  the  steady-state  sample  covariance 
of  e^  is  the  estimate  of  t.  Since  x*  for  this  approach  starts  from 
zero,  the  number  of  iterations  to  get  the  estimate  of  x  depends  on  the 
value  of  x.  If  x  is  close  to  zero  or  a  small  value,  then  x*  can  be 
estimated  in  a  few  iterations.  The  maximum  number  of  iterations  of 
this  approach  Is  equal  to  NT  (when  x  is  equal  to  T).  This  maximum  number 
of  Iterations  Is  greater  than  the  maximum  number  of  Iterations  of  the 
previous  approach.  For  example,  if  NT  is  equal  to  2^,  then  the  maximum 
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number  of  iterations  of  this  approach  is  equal  to  2V  while  that  of 

the  previous  approach  is  only  equal  to  V. 

As  an  example,  consider  the  system  in  Figure  2  of  Section  II 

with  a  second-order  plant  and  a  first-order  controller.  The  external 

input  Wp  is  Gaussian  white  noise,  which  has  zero  mean  and  variance  =  1. 

Let's  assume  t  to  be  one  of  the  values  0,  _J_,  .  .  .  ,  T  and  let  NT 

NT 

be  50.  A  FORTRAN  program  to  simulate  the  entire  closed-loop  system 

in  Figure  7  and  to  implement  the  algorithm  for  estimating  x  of  the 

above  example  is  in  Appendix  F.  All  arrays  of  this  program  are  the 

same  as  those  of  the  program  of  the  new  model. 

The  system  in  Figure  7  is  simulated  by  the  software  in  Appendix 

F.  The  software  which  implements  the  algorithm  for  estimating  t, 

will  wait  until  the  state  observers  (x*)  equal  the  state  variables 

(xp)  and  this  system  is  in  steady  state.  From  Appendix  A  (this  system 

is  in  steady  state  at  the  time  200T)  and  Appendix  C  (the  state  observers 

equal  the  state  variables  at  time  approximate  equals  40T),  the  software 

will  wait  for  200T,  then  this  software  starts  to  implement  the  algorithm 

for  estimate  x.  As  in  the  example  in  Section  II  and  Section  III,  N 

is  given  the  value  100  in  the  equation  for  calculating  the  sample 

covariance  of  the  errors.  With  t  which  is  supposed  to  be  equal  to 

3T  ,  this  software  computes  the  estimate  of  this  t  which  is  equal  to 
NT 

XL  and  the  number  of  iterations  which  is  equal  to  5. 

NT 

In  this  example,  t  Is  assumed  to  be  one  of  the  values  0,  _T, 

NT 
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2T,  .  .  .  ,  T,  and  x*  is  restricted  to  be  any  value  among  0,  J[, 

NT  NT 

.  .  .  ,  T.  Therefore,  t  can  be  estimated  exactly.  In  general, 

xTwill  be  between  0  and  T,  but  it  may  not  be  one  of  the  above  values. 

Thus,  the  estimate  of  x  will  generally  not  equal  to  x  and  the  maximum 

difference  between  x  and  x*  is  T.  The  numerical  values  of  this 

NT 

difference  can  be  reduced  by  increasing  NT. 

Since  the  time  skew  varies  with  time,  it  should  be  estimated 
at  regular,  short  interval  of  times.  Thus,  the  execution  time  of  the 
algorithm  in  this  subsection  for  estimating  x  should  be  very  small. 

From  the  above  example,  the  total  execution  time  for  estimating 
x  is  the  sum  of  the  time  for  computing  Pe/\ss  and  the  time  for  executing 
the  algorithm  for  estimating  x.  The  time  for  executing  the  algorithm 
for  estimating  x  can  not  be  reduced  but  the  time  for  computing  PfiAss 
can  be  reduced  only  by  decreasing  the  number  of  values  of  e^.  (The 
difference  between  y£2(tk  +  t*)  and  yC2(t|<  +  x).) 

Since  the  sample  covariance  of  e^  in  equation  3-19  is  the  estimated 
value  of  the  actual  covariance  of  eA,  the  difference  between  the 
estimated  value  and  the  actual  value  depends  on  the  number  of  values 
of  eA  used.  By  the  'Law  of  Large  Number'  of  Probability  and  Statistics, 
if  N,  which  is  the  number  of  random  variables  of  e^  used,  is  large, 
there  is  a  high  probability  that  e^  (the  sample  mean  of  e^)  will  be 
closed  to  the  actual  mean  of  these  random  variables.  Then,  there  is 
high  probability  that  the  sample  covariance  of  e^  will  be  closed  to 
the  actual  mean,  too.  Otherwise,  if  N  is  small,  the  estimated  value 
may  diverge  from  the  actual  value.  According  to  the  algorithm  for 
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estimating  t,  the  accuracy  of  this  algorithm  only  depends  on  how 
smooth  the  curve  PeAss  is.  Then  the  covergence  or  divergence  of  the 
estimated  value  (the  sample  covariance  of  eA)  to  the  actual  value 
(the  covariance  of  e^)  does  not  concern  to  this  algorithm.  The  value 
of  N  can  be  selected  as  small  as  the  curve  of  PeAss  is  plotted 

as  a  function  of  t  and  t*  is  still  smooth.  However,  with  N  =  10,  the 
curve  of  PeA$$  is  as  smooth  as  the  curve  of  PeAs$  w1th  N  =  *00* 

2.  CHARACTERISTICS  OF  THE  TIME  SKEW 

By  assumption,  the  time  skew  is  constant  over  a  short  period  of 
time.  However,  in  reality  the  value  of  time  skew  will  change  slowly 
with  time.  As  discussed  in  Reference  1,  the  time  skew  can  be  assumed 
to  vary  linearity  with  time  and  it  is  equal  to  T  at  the  time  which  the 
two  samplers  return  to  synchronism.  This  characteristic  variation  of 
time  skew  is  shown  in  Figure  12. 

In  the  algorithm  for  estimating  t  in  last  subsection,  t  is 
assumed  to  be  changing  very  slowly.  However,  this  algorithm  cannot 
estimate  t  at  that  time  at  which  t  changes  from  T  to  zero.  Another 
means  for  estimating  t  is  needed  for  this  time  interval. 

After  the  first  estimate  of  t  is  obtained,  the  characteristic 

curve  in  Figure  13  can  be  drawn.  Let  0A  in  Figure  13  represent  the 

time  that  corresponds  to  the  first  t*,  which  is  represented  by  AB. 

Since  the  time  skew  varies  linearly  with  time,  the  estimate  of  the 

time  skew  also  varies  linearly  with  time.  The  characteristic  of  x * 

can  be  drawn  by  using  the  slope  M  and  the  time  which  x*  equals  T 

W 

Is  the  approximate  value  which  the  two  sample  periods  coincide.  Thus, 
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an  approximation  to  the  time  at  which  two  sample  periods  coincide 

(point  a  in  Fiqure  13)  is  qiven  by  T  x  AB  .  Thus,  after  the  first 

OB 

t*  has  been  estimated,  the  time  at  which  two  sample  periods  coincide 
is  estimated  from  this  latter  equation. 
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SECTION  V 


OPERATION  OF  THE  NEW  MODEL 


This  section  describes  the  asynchronous  operation  of  a  digital 
flight  control  system  using  the  technique  in  Section  III.  The  output 
of  the  first  controller  in  the  model  used  is  always  the  input  to  the 
plant.  The  outputs  of  both  controllers  are  sent  to  the  monitor,  which 
will  compare  these  signals  using  the  tolerance  value  from  the  technique 
in  Section  IV. 

1 .  DUAL-REDUNDANT  DIGITAL  FLIGHT  CONTROL  SYSTEM. 

A  model  for  a  dual -redundant  digital  flight  control  system  is 
shown  in  Figure  14.  The  output  of  the  plant  is  sampled  by  each  of  the 
controllers,  using  a  common  sample  period  but  having  a  fixed  time  skew 
between  them.  The  output  of  the  first  controller  serves  as  the  input 
to  the  plant.  The  outputs  of  the  controllers  go  to  the  monitor  and  the 
monitor  will  first  isolate  the  bad  signal  and  then  select  or  calculate 
the  best  signal  from  the  remaining  good  signals. 

Figure  15  shows  the  details  of  the  monitor.  All  the  previous 
values  of  output  of  the  plant  (yp(tk_-| ) ) ,  the  output  of  channel  1 
(yci (^k-l ) ) *  an<*  externa^  input  to  the  plant  (wp( t^.-j ) )  go  to  the 

block  named  OBSERVER,  which  is  a  subsystem  designed  to  estimate  the 
state  variables  of  the  plant.  The  observer  produces  XpUfc),  an  estimate 
of  XpUfc).  xpU|<  +  T*)«  an  estimate  of  xp(t|C  +  t),  is  computed  from 
xp(*k)»  ycl(V»  and  T*‘  Then  yc2(tfc  +  t*),  an  estimate  of  yc2(tk  ♦  t*) 
that  is  used  to  calculate  the  Inherent  error,  is  computed  from 


FIGURE  14  :  VOTING  OF  DUAL- REDUNDANT  CHANNEL 


X*(tk  +  T*)  and  the  state  variables  of  the  2nd  DIGITAL  CONTROLLER  #1. 

P  K 

After  the  system  in  Figure  14  estimates  x,  both  +  T*) 

and  yc2(tk  +  t)  are  sent  to  the  monitor  logic.  If  the  sample  covariance 
of  the  difference  between  y£2(tk  +  T*)  anc*  y^C^k  +  ’s  9reater 
than  the  tolerance  value  which  is  discussed  in  Section  IV,  then 
channel  failures  are  probably  present.  The  monitor  logic  will  isolate 
the  bad  signal  (yc-j  or  y^)  from  the  channel  failure. 

As  discussed  at  the  end  of  Section  IV,  the  execution  time  for 
estimating  x  can  be  reduced  by  decreasing  N  in  the  equation  of  the 
sample  covariance  of  e^  (eq.  3-19).  Section  IV  also  shows  that  there 
is  no  difference  in  the  estimate  of  x  for  N  equal  to  100  or  10.  Thus, 

N  in  this  section  is  selected  to  be  equal  to  10  for  the  purpose  of 
reducing  the  execution  time  above. 

To  estimate  x,  the  2nd  DIGITAL  CONTROLLER  #1  executing  the  algorithm 
in  Section  IV  for  estimating  x  will  wait  50T  seconds  (T  =  0.0125  second) 
until  the  state  observer  (xp)  is  equal  to  the  state  variables  (xp) 

(The  details  are  in  Appendix  C.).  Then  the  2nd  DIGITAL  CONTROLLER  #1 
will  store  the  next  ten  values  of  e^  for  estimating  x.  After  50T 
seconds,  x  can  be  estimated  in  the  time  required  for  obtaining  10 
values  of  e^  plus  the  time  for  estimating  x.  Since  the  ten  values  of 
e^  (from  51 T  second  to  60T  second)  are  used  to  estimate  x,  then  this 
x*  is  the  estimated  value  of  x  between  51T  second  and  60T  second. 

While  the  2nd  DIGITAL  CONTROLLER  #1  is  In  the  process  of  estimating 
x,  ycl  and  y ^  are  always  sent  to  the  monitor.  Then,  for  the  time 
interval  o  to  the  time  at  which  the  first  x*  is  estimated  (50T  seconds 
are  required  for  the  state  observers  to  equal  to  the  state  variables.). 
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the  monitor  logic  should  be  disabled. 

As  in  Reference  1,  let  P  be  equal  to  the  period  required  for 
T,  and  to  return  to  synchronism  or 


Let  e  be  the  fractional  error  between  the  clock  crystals  controlling 
the  separate  processors;  i.e.. 


h 

then  P  =  T2  =  T1  ”^1  ^1 

e  e  e 

With  Tj  =  0.0125  second  and  e  =  0.017,  it  requires  10,000  samples 
of  7}  and  T2  to  return  to  synchronism.  Thus,  t  changes  very  slowly 
with  time,  then  the  t*  in  the  previous  time  interval  can  reasonably 
be  assumed  to  be  the  estimated  time  skew  over  the  current  time  interval. 

As  discussed  in  Section  IV,  the  time  at  which  the  two  sample 
periods  return  to  synchronism  is  important  because  t  at  this  time 
changes  abruptly  from  T  to  0.  For  the  time  at  which  x  is  equal  to  T 
or  0,  it  can  be  estimated  by  the  technique  described  in  the  last  sub¬ 
section  of  Section  IV. 

After  t*  is  estimated,  the  monitor  will  compare  y^Ct^  +  T*) 
and  yc2(t^  +  T)*  If  the  difference  of  these  two  signals  is  greater 
than  the  given  tolerance  value  discussed  in  Section  IV,  then  a  channel 
failure  has  occurred. 
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2.  TRI-REDUNDANT  DIGITAL  FLIGHT  CONTROL  SYSTEM 

A  model  for  a  tri -redundant  digital  flight  control  system  is 
shown  in  Figure  16.  The  output  of  the  plant  is  sampled  by  each  of 
the  controllers,  using  a  common  sample  period  but  having  two  fixed 
time  skews  x^  and  x2  between  channel  1  and  2,  and  channel  1  and  3, 
respectively.  The  output  of  the  first  controller  serves  as  the  input 
to  the  plant  and  the  outputs  of  these  three  controllers  go  to  the 
monitor,  as  in  the  model  of  the  previous  subsection. 

Figure  17  shows  the  details  of  the  monitor.  The  output  of  the 
first  controller  is  used  to  calculate  y^ft^  +  x*),  an  estimate  of 
yC2(tk  +  Tjh  y^k  +  T2) *  an  estimate  of  yc3(tk  +  t2),  where 
t*  and  t*  are  the  estimated  values  of  x^  and  x£,  respectively.  Using 
the  same  technique  as  discussed  in  the  last  subsection,  x-j  can  be 
estimated  from  y£3(tk  +  *2)  and  +  T2^-  since  the  maximum 

differences  between  x*  and  x-j  ,  and  x£  and  x ^  are  equal,  then  the  tol¬ 
erance  values  between  channel  1  and  2,  and  channel  1  and  3  are  equal. 

In  comparing  channel  2  and  3,  the  tolerance  value  is  approximately  equal 
to  the  maximum  covariance  of  the  difference  between  +  Tl)  and 

y*3(tk  +  X2>)-  Thus,  if  the  covariance  of  the  difference  between 
yc2(tk  +  x-j )  and  yc3Uk  +  *2)  1s  greater  than  the  tolerance  value 
between  channel  2  and  3  described  above,  the  malfunction  must  have 
occurred  in  channel  2  or  channel  3. 

The  model  for  more  than  three  channels  can  be  described  in  a 
manner  similar  to  the  model  for  the  three-channel  system.  For  example, 
consider  a  four-channel  system.  This  system  has  three  time  skews: 
the  time  skew  between  channels  1  and  2  (x^),  the  time  skew  between 
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FIGURE  16  :  VOTING  OF  TRI- REDUNDANT  CHANNEL 
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yci(tK-i 


channels  1  and  3  (x^),  and  the  time  skew  between  channels  1  and  4 
(t^).  The  2nd  DIGITAL  CONTROLLER  will  estimate  these  three  time 
skews  for  the  estimated  value  of  the  output  of  channel  2,  the  estimated 
value  of  the  output  of  channel  3,  and  the  estimated  value  of  the  out¬ 
put  of  channel  4.  As  in  the  three-channel  system,  the  maximum  differ¬ 
ences  between  x^  and  x*.  the  estimate  of  x^ ,  and  T2»  ^e  estimate 
of  X2»  X3  and  x^,  the  estimate  of  x^  are  equal  to  T/NT.  Then  the 
tolerance  values  of  channels  1  and  2,  channel  1  and  3,  and  channel  1 
and  4  are  equal  to  the  maximum  sample  covariance  between  the  estimated 
value  and  the  actual  value  of  any  one  of  the  outputs  of  the  controller 
(the  maximum  sample  covariance  of  the  differences  between  the  estimated 
and  the  actual  values  of  the  output  of  the  second  channel,  the  estimated 
and  the  actual  values  of  the  output  of  the  third  channel  and  the  estimated 
and  the  actual  values  of  the  fourth  channel  are  equal).  For  channels  2 
and  3,  the  tolerance  value  is  equal  to  the  maximum  covariance  of  the 
difference  between  the  estimated  values  of  channel  2  and  3.  Similarly, 
the  tolerance  value  of  channels  2  and  4,  and  the  tolerance  value  of 
channel  2  and  4  are  equal  to  the  maximum  covariance  of  the  differences 
between  the  estimated  values  of  channel  2  and  4  and  channel  3  and  4, 
respectively. 
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SECTION  VI 


COMPARISON  BETWEEN  BASIC  MODEL 
AND  NEW  MODEL 

The  basic  model  can  distinguish  inherent  errors  from  the  errors 

induced  by  channel  failures  by  using  the  maximum  steady-state  sample 

covariance  of  e^  as  the  tolerance  value.  If  the  measured  covariance 

of  e.  is  greater  than  this  tolerance  value,  then  the  monitor  indicates 
A 

a  channel  failure.  Otherwise,  if  the  measured  covariance  of  e^  is 
less  than  this  tolerance  value,  the  monitor  indicates  that  only  inherent 
errors  are  present.  However,  this  tolerance  value  is  not  the  best 
value  to  use  to  distinguish  inherent  errors  from  errors  induced  by 
a  channel  failure.  If  the  measured  covariance  of  e^  is  greater 
than  the  covariance  of  e^  of  the  present  t  but  is  less  than  the  maximum 
steady-state  sample  covariance  of  eft  (t  *  T),  then  the  basic  model  would 
not  indicate  the  channel  failure.  To  reduce  this  deficiency,  the  new 
model  computes  a  tolerance  value  equal  to  the  maximum  steady-state 
sample  covariance  of  the  difference  between  y^(tk  +  T*)  and  +  *)• 

Since  this  tolerance  value  is  very  small  when  it  is  compared  with  the 
tolerance  value  of  the  basic  model,  the  deficiency  of  the  basic  model 
discussed  above  can  be  reduced. 

The  number  of  tolerance  values  of  the  new  model  depends  on  the 
number  of  channels.  There  is  one  tolerance  value  for  two  channels, 
two  tolerance  values  for  three  channels,  four  tolerance  values  for 
four  channels,  and  so  on.  For  the  basic  model,  there  is  only  one 
tolerance  value  for  any  number  of  channels. 
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The  next  comparison  is  the  hardware  structure.  Both  models  (the 
basic  model  and  the  new  model)  require  a  computer  to  calculate  the 
covariance  of  e^..  But  the  software  of  the  new  model  is  more  complicated 
than  that  of  the  basic  model.  The  software  of  the  new  model  is  used  to 
compute  x*,  an  estimate  of  xp,  Xp(tk  +  r*),  an  estimate  of  Xp(tk  +  t), 
y^Uk  +  T*)»  an  estimate  of  yc2(tk  +  O.  t*,  an  estimate  of  t,  and 
the  sample  covariance  of  (The  difference  between  ycl(tk)  and 
yc2(tk  +  t)  is  eA  of  the  basic  model.  The  difference  between  *  '*> 

and  yC2(t|<  +  T)  1S  eA  the  new  model.). 
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SECTION  VII 

SUMMARY  AND  RECOMMENDATIONS 

Two  models  for  the  asynchronous  digital  flight  control  system 
are  described  in  this  report.  The  basic  model  is  almost  the  same  as 
the  basic  model  in  Reference  1  and  2  except  that  the  input  to  the 
plant  is  the  difference  between  the  external  input  and  the  output  of 
the  controller.  However,  the  input  to  the  plant  is  always  the  output 
of  the  first  controller. 

As  discussed  in  Reference  1,  this  basic  model  is  roughly  equivalent 
to  the  voter  named  'median  select*  (the  upper  median  for  a  four- 
channel  system  or  the  lower  median  for  a  three-  or  four-channel  system 
is  used  as  the  voter  output)  when  the  channel  outputs  are  either 
monotonically  increasing  or  decreasing  in  time.  Figure  18  illustrates 
the  example  of  an  asynchronous,  dual -redundant  digital  flight  control 
system  which  produces  a  monotonically  increasing  output.  Channel  1 
produces  the  sampled  outputs  at  times  t^,  t^  ,  .  .  .  ,  for  k  =  0, 

1,  .  .  .  and  channel  2  produces  the  sampled  outputs  at  times  t^, 

tk  +  t,  tk+j  +  t,  .  .  .  ,  for  k  =  0,  1 . A  comparison  monitor 

which  compares  the  magnitudes  of  the  outputs  of  the  two  channels  will 
observe  differences  illustrated  as  e^  and  eg  in  Figure  18. 

Since  the  pilot's  command  or  wind  -  gust  which  changes  all  the 
time  is  the  external  input,  then  the  random  signal  (Gaussian  white 
noise)  is  chosen  to  be  the  external  input  of  the  models.  A  comparison 
monitor  in  this  case  (Gaussian  white  noise  is  the  external  input.)  will 
compare  the  signal  by  using  the  covariance  of  e^  instead  of  e^. 
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FIGURE  18  INHERENT  ERRORS  IN  ASYNCHRONOUS  OPERATION 


According  to  the  results  of  the  basic  model  in  Section  II,  the 
steady-state  sample  covariance  of  e^  (PeAs$)  is  largest  when  the  times 
at  which  ycj  (output  of  the  first  controller)  and  yC2  (output  of  the 
second  controller)  change  are  farthest  apart.  Then  the  tolerance 
value  of  this  model  is  equal  to  the  maximum  steady-state  sample 
covariance  of  the  difference  between  y^Ct^)  and  y  +  0  (when 
T  =  T). 

The  new  model  described  in  Section  III  is  an  extension  of  the 
basic  model  in  Section  II.  As  in  the  basic  model,  the  external  input 
to  this  model  is  a  Gaussian  white  noise.  This  model  tries  to  reduce 
the  deficiency  of  the  basic  model  described  at  the  beginning  of  Section 
III  by  decreasing  the  tolerance  value.  There  are  two  functions  per¬ 
formed  by  the  first  channel  of  this  model.  The  first  function  is  to 
compute  the  control  output  to  the  plant.  The  second  function  is  to 
compute  a  signal  used  to  calculate  the  inherent  error;  this  signal  is 
an  estimated  value  of  the  output  of  the  second  channel. 

According  to  the  results  of  the  new  model  in  Section  III,  the 
steady-state  sample  covariance  of  the  difference  between  y^U^  +  t*) 
and  y^tfc  +  t)  is  directly  proportional  to  the  difference  between 
t  and  t*.  If  t*  equals  t,  the  difference  between  y^ft^  +  T*)  and 
yc2(t|c  ♦  t)  Is  equal  to  zero.  Then  the  tolerance  value  of  the  new 
model  depends  on  the  difference  between  t  and  t*. 

The  algorithm  in  Section  IV  for  estimating  t  Is  based  on  the 

results  above.  Let  t  and  t*  be  one  of  the  values  0,  .  .  .  ,  T. 

NT 

Then  this  algorithm  computes  the  steady-state  sample  covariance  of  the 
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difference  between  y* g(tk  +  x*)  and  y^t^  +  T)  when  T*  is  the  midvalue 

between  the  first  lower  limit  (0)  and  the  first  upper  limit  (T).  If 

the  covariance  of  the  difference  between  y£2(t|c  +  T*)  and  yc2(t|c  + 

of  this  first  x*  is  less  than  that  of  the  next  value  of  this  first 

x*.  then  must  be  between  0  and  the  new  upper  limit  (the  midvalue). 

Otherwise,  if  that  of  this  first  x*  is  greater  than  that  of  the  next 

value  of  this  first  x*,  then  x  must  be  between  the  new  lower  limit  (the 

midvalue)  and  T.  By  using  this  scheme,  the  interval  containing  x  can 

be  reduced  by  half  for  each  iteration.  The  algorithm  will  repeat  this 

technique  until  the  current  midvalue  differs  from  the  previous  midvalue 

by  _T.  In  the  last  interval,  any  values  of  x*  which  corresponds  to  the 
NT 

smallest  steady-state  sample  covariance  of  the  difference  between 

y*2(tk  +  x*)  and  yc2(tk  +  x)  is  the  estimate  of  x.  In  general,  x  may 

not  be  one  of  these  values  0,  _J,  .  .  .  ,  T,  then  the  maximum  difference 

NT 

between  x  and  x*  is  _J\ 

NT 

The  results  from  Section  III  and  Section  IV  can  be  applied  to 
the  asynchronous  operation  of  digital  flight  control  system.  For  a 
two-channel  system,  the  DIGITAL  CONTROLLER  #1  in  the  monitor  will  wait 
until  the  state  observers  (x*)  equal  the  state  variables  (xp).  Then 
the  DIGITAL  CONTROLLER  #1  will  execute  the  algorithm  for  estimating 
x  by  using  the  next  10  values  of  the  difference  between  y$2(tk  +  x*) 
and  yc2(tk  +  x).  During  the  time  DIGITAL  CONTROLLER  #1  is  estimating  the 
first  x*,  the  monitor  will  not  compare  the  signals  ^c2^k  +  T*)  and 
yc2 (tk  +  T))*  After  the  first  x*  is  estimated,  the  monitor  will  compare 
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the  signal  by  using  the  first  t*  until  the  new  t*  is  estimated.  Then 
the  monitor  will  compare  the  signals  by  using  the  new  t*  and  so  on. 
However,  during  the  time  in  which  the  first  t*  is  estimating  (after  the 
state  observers  (xp)  equal  the  state  variables  (xp)),  the  monitor  can 
compare  the  signals  by  using  the  tolerance  value  which  is  equal  to  the 
sample  covariance  when  t*  =  0  and  t  =  T.  t*  can  be  estimated  for  every 
time  period  except  the  time  at  which  t  changes  from  its  maximum  value 
to  its  minimum  value  (The  details  are  in  Section  IV.). 

For  a  system  with  three  or  more  channels,  the  number  of  time  skews 
depends  on  the  number  of  channels.  For  example,  there  are  two  time 
skews  for  the  three-channel  system;  r^,  the  time  skew  between  channels 

1  and  2;  and  t£,  the  time  skew  between  channels  1  and  3.  For  the 
four-channel  system,  there  is  one  more  time  skew:  t3,  the  time  skew 
between  channels  1  and  4.  All  these  time  skews  can  be  estimated  by 
using  the  same  technique  as  in  the  two-channel  system.  After  these 
time  skews  are  estimated,  the  monitor  will  compare  channels  1  and  2, 
channels  1  and  3,  and  so  on  by  using  the  estimated  output  of  channel 

2  and  the  actual  output  of  channel  2  for  channels  1  and  2  and  so  on.  If 
the  sample  covariance  of  the  difference  between  the  estimated  and  the 
actual  values  of  channels  1  and  2  or  channels  1  and  3  and  so  on  is  greater 
than  the  tolerance  value  of  the  new  model,  then  the  channel  failure  is 
occurred.  For  a  pair  of  the  channels  2,  3,  ...  ,  the  monitor  will 
compare  a  pair  of  the  channels  by  using  the  estimated  outputs  of  that 
pair  of  the  channels.  For  example,  the  monitor  will  compare  channels 

2  and  3  by  using  the  tolerance  value  which  is  the  maximum  covariance 
of  the  difference  between  the  estimated  values  of  channels  2  and  3.  If 
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the  sample  covariance  of  the  actual  outputs  of  channels  2  and  3  is 
greater  than  the  above  tolerance  value,  then  a  channel  failure  has 
occurred. 

The  last  section  describes  the  comparison  between  the  basic 
model  and  the  new  model.  The  disadvantage  of  the  basic  model  is  that 
the  basic  model  would  not  indicate  the  channel  failure  although  the 
sample  covariance  is  greater  than  the  sample  covariance  of  the  present 
t.  The  new  model  can  reduce  this  deficiency  by  reducing  the  tolerance 
value.  However,  the  new  model  requires  more  hardware  and  is  more 
complicated  to  simulate. 

It  is  recommended  that  the  work  be  continued  to  accomplish  the 
following: 

1.  The  software  in  this  report  can  only  implement  the  example 
in  Figure  2.  Then  the  software  for  simulating  the  models  should  be 
developed  to  be  the  general  software.  Thus,  software  for  simulating 
a  class  of  systems  is  needed. 

2.  Increase  the  complexity  and  generality  of  the  models  and 
covariance  analysis  to  include  such  features  as  multirate  sampling, 
computational  delays,  processor  word-length  effects,  sensor  noise  and 
additional  voter  algorithms. 

3.  To  reduce  the  complexity  of  the  asynchronous  operation  of  a 
digital  flight  control  system  using  the  technique  in  Section  IV,  a 
model  of  a  new  algorithm  for  a  dual -redundant  system  is  shown  in  Figure 
19.  In  this  model,  t*  is  constant  and  equal  to  T/2.  After  the  state 
observers  are  equal  to  the  state  variables,  the  signals  4<*k  ♦ 1/2 1- 
and  yC2(*k  +  T)  *re  sent  t0  the  won^or  logic.  The  tolerance  value  to 
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use  in  this  case  is  equal  to  the  sample  covariance  of  the  difference 
between  y^U^  +  T/2)  and  yC2 ( +  0).  This  algorithm  is  simpler  than 
the  algorithm  described  in  Section  V  because  t*  is  constant.  But  the 
tolerance  value  of  this  new  algorithm  is  larger  than  the  tolerance 
value  of  the  algorithm  in  Section  V.  However,  the  tolerance  value  of 
this  algorithm  is  less  than  the  tolerance  value  of  the  basic  model  by 
half. 

A  model  of  a  new  algorithm  for  a  tri -redundant  system  is  shown  in 
Figure  20.  In  this  system,  the  tolerance  value  between  channels  1  and 
2,  and  channels  1  and  3  are  equal  to  the  sample  covariance  of  the 
difference  between  y^Uk  +  T/2)  and  y^t^  +  T)-  It  is  not  necessary 
to  estimate  yc3(tk  +  T/2)  because  yC3(tk  +  T/2)  is  equal  to  y^t^  +  T/2). 
However,  the  tolerance  value  of  channels  2  and  3  is  the  same  as  in  the 
basic  model. 

4.  The  number  of  values  of  e^  for  computing  the  sample  covariance 
should  be  studied  with  the  objective  of  using  as  few  values  as  needed 
for  a  reasonable  reduction  in  the  variability  of  the  estimate. 
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APPENDIX  A 

SOFTWARE  FOR  THE  BASIC  MODEL 

1.  FLOWCHART  ANO  DESCRIPTION  OF  MAJOR  COMPONENTS  AND  SUBROUTINE 

As  in  Reference  1,  the  main  program  of  the  software  for  the 
basic  model  is*  called  PROGRAM  SKEW.  Its  major  computational  tasks 
are  to  develop  the  state  variable  model  of  the  complete  closed-loop 
system  and  to  compute  the  controller  outputs,  the  errors  and  eg, 
and  the  steady-state  covariances  of  the  states. 

A  flowchart  for  the  program  appears  in  Figure  21.  The  blocks 
in  this  figure  correspond  to  the  clearly  identified  components  of  the 
main  program. 

The  first  block  shows  the  data  input.  The  variables  are  self- 
explanatory  except  for  the  quantities  NT,  NTAU,  and  NT2.  Since  a 
numerical  integration  is  required  to  compute  tp(x) ,  and  <|>(T),  it  is 
necessary  to  quantize  the  time  interval  [0,T].  The  user  specifies 
the  degree  of  quantization  by  specifying  NT,  the  number  of  subintervals 
in  [0,T]  which  are  to  be  used  in  the  computation.  For  convenience, 
the  subintervals  of  [0,T]  are  designated  1,2,  .  .  .  to  NT1;  where  NT1 
is  equal  to  NT  +  1.  NTAU  represents  the  value  of  which  Is  computed 
within  the  program  as 

ri  =  (NTA^  -  1)  *  T  B-l 

(NT1  -  1) 
and 

NTAU,  *  ITAU  *  NT  +  1  B-2 

1  NT? 
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START 


Read  in  heading  and  description  of  the  plant  and 
the  controller ;%p."p.nwp,nop,nc,Ap,Blp,Cp,Fc,Gc, 


HC,T,  NT,  and  NT2 


Compute  the  Gaussian  white  noise  by  using 
SUBROUTINE  RANDU 


Calculate  and  *(T)  by  using  eq.2-4 


Calculate  ¥ (t  j) ,  and  ¥ (T)  by  using  the  trapezoidal 
rule  for  numerical  integration  _ 


Calculate  xc2^tK+r^'  yc2^tK-l+r^  using  eqs.  2-10 
and  2-11 


Calculate  xci(tR+1),  yci(tR)  by  using  eqs.  2-8 
and  2-9  _ 


Calculate  P^gg.  and  PeBs8  using  eqs.  2-23  and  2-24 


'r'  *--•■*»  -•■  ■■■■ 


for  ITAU  =  0,  1,  2,  .  .  .  ,  NT2 

where  NT_  must  be  an  integer. 

NT2 

In  block  two,  the  gaussian  noise  is  computed  by  using  the  sub¬ 
routine  named  RANOU.  This  subroutine  is  available  in  most  IBM-based 
computer  systems. 

The  third  block  specifies  the  calculation  of  ♦(r<)  and  ♦(T). 

The  computations  require  4  from  block  2,  so  that  the  required  numerical 
integrations  can  be  performed.  The  numerical  integrations  use  the 
trapezoidal  approximation. 

In  block  five  of  Figure  21,  the  second  controller  state  variable 
and  the  second  controller  yc2  are  calculated  by  using  equations 

xc2^k  +  =  Fcxc2^k-1  +  +  Gcuc2^k-1  + 

yC2^k-l  +  T)  *  Hcxc2{tk-1  +  T>  +  EcUc2(tk-l  +  t) 

In  block  six,  the  first  controller  state  variable  xc^  and  the 
first  controller  output  y£^  are  calculated  by  using  equations 

xcl^tk+l  ^  =  Fcxcl^k^  +  Gcucl(V 

ycl(V  -  WV  +  EcUcA> 

In  block  seven,  the  inherent  errors  eA  and  eg  are  calculated 
by  using  equations 

eA(t)  *  ycl(V  '  yc2(tk  +  »> 


I 
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r 


fbrtk  +  T<t<tk+1,0<T<T.  k-O.  1.  .  .  . 

®B(t)  *  ycl(tk+l*  -  yc2(tk  +  T> 

for  *k+l  -  1  *  Vl  *  '•  0  <  T  -  T-  k  ’  °*  ’•  ’  •  • 

The  final  set  of  computations  is  given  in  block  eight.  The  steady 
state  covariance  PEASS  and  PEBSS  are  calculated  by  using  the  equations 

peAss  =  jjf  ^  teAi  "  FA^CAi  *  FA^ 

PeBss  =  J  iJ1  CeBi  ’  eB^eBi  "  V 

where 

e^  and  e^.:  are  the  errors  when  the  system  is  in  steady 
state 
and 

eA  and  are  the  sample  means  of  N  samples  of  e ^  and  egi 

2.  INSTRUCTIONS  FOR  USING  THE  PROGRAM 

The  first  data  card  is  used  to  provide  a  message  which  will  be 
printed  at  the  top  of  a  new  page  of  output.  The  next  card  specifies 
NP,  NUP,  NWP,  NOP,  and  NC  using  the  format  (513).  These  quantities 
are  the  actual  dimensions  of  the  plant  and  controller.  Next,  the 
matrices  Ap,  Bp,  Cp,  Ec>  Fc»  Gc,  and  Hc  are  specified  In  succession, 
one  row  and  one  card  at  a  time,  using  the  FORMAT  (F10.4,  213).  The 
next  card  specifies  T,  NT,  and  NT2  using  FORMAT  (F10.4,  315).  Table 
1  shows  what  values  to  assign  to  the  given  arrays. 

The  computer  program  listing  for  PROGRAM  SKEW  of  the  basic  model 
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■>  Nr»  f-.TT..- 


written  in  FORTRAN  appears  in  appendix  B.  PEASS  and  PEBSS  for  each 
value  of  t  are  shown  at  the  end  of  the  listing. 


TABLE  1 

REQUIRED  DIMENSIONS  OF  ALL  ARRAYS 


AP  (NP,  NP) 

YP  (NOP) 

BP  (NP,  NWP) 

XPTAU  (NP) 

CP  (NOP,  NP) 

YPTAU  (NOP) 

FC  (NC,  NC) 

UP  (NUP) 

GC  (NC,  NOP) 

W  (1000) 

HC  (NUP,  NC) 

YC1  (NUP) 

EC  (NUP,  NOP) 

XC1  (NC) 

ECCP  (NUP,  NP) 

YC2  (NC) 

PHST1  (NP,  NP) 

XC2  (NC) 

PHTAU  (NP,  NP) 

El  (NUP) 

PHTAU1  (NP,  NP) 

E2  (NUP) 

PS ST  (NP,  NUP) 

EA  (1000) 

PSTAU  (NP,  NUP) 

EB  (1000) 

PSTAU1  (NP,  NUP) 

PEASS  (30) 

XP  (NP) 

PEBSS  (30) 

Since  the  example  is  in  the  steady-state  at  time  approximately  equal 
200T  (T  =  0.0125  second),  then  the  first  value  of  i  in  equation  (2-23), 
(2-24),  (2-25),  and  (2-26)  is  201  and  let's  assume  the  value  of  N  in 
these  equations  equals  100. 


APPENDIX  B 

COMPUTER  PROGRAM  LISTING 
FOR 

PROGRAM  SKEW  AND  EXAMPLE 
OF  OUTPUT  WRITTEN  IN  FORTRAN 
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COMMON  EL.EMX  *  MAXI  *  MAX  J 

DIMENSION  AF'  (2  f  2 )  f  B1P< 2 f  2 )  *  CF‘<  1 » 2 )  t  FC  <  2 f 2 )  »GC  (  2f  1 )  fPHIT  <  4 f 4 ) f 
1  HC(2f2)fEC(2f1)  iPHIT1(4f4f  101 )  » PSIT1  <  4  ,  4 )  ,PHTAU<  4f  4 >  fPSTAU<4f4) f 

4  INEiEX  (4 )  f W<  4000  )  f  F’S  <  4  f  4  >  *YC1<2)  »F‘SIT  <4f4)fXW3<2)f 

5  YC2  (2  >  fEI  <  2)  f E2< 2 )  f XF‘< 2 )  f  XPTAU<  2 )fXP1<2)fYP<2)f YPTAUl 2 ) f 

6  XC1(2) f XC2<2) f AM<4f 4) r PTC 4f 4 ) »P1 ( 4f4 ) , D1 <4 » 4) » 

7  D2(4f4) fD3<  4) fXWI ( 2)  f XW2< 2) »ECCP<4*4) »D<4»4> » 

8  PEASS  ( 50 )  f PEBSS  <  50 )  fEA(IIOO) fEB< 1100) 

CCCCC  PROVIDE  MAXIMA  FOR  CALLED  ARRAYS 

NPM  "  2 
NUF'M  -  2 
NWPM  =  2 
NOPM  =  1 
NCM  =  2 

NHM  =  NPM  +  NUPM 
NFM  «  NPM  +  2*NCM 
NRRM  -  2#NPM  +4 

C 

C  READ  INPUT  DATA 

C  ###*****#**#*###*##***#*##**##*###*###**#*********#******#******* 

C 

WRITE <6f 899) 

899  FORMAT ( ' 1' ) 

100  READ(5r 900)  ID 
900  FORMAT (20A4) 

WRITE (6 r 902)  ID 
902  FORMAT ( '1 ' f20A4) 

READ  <  5 f 90 6 ) NP f NUP f NWP  f  NOP  f  NC 
906  FORMAT (513) 

WRITE ( 6 f 908 )  NP f NUP r NWP f NOP f NC 
908  FORMAT ('ONG.  OF  PLANT  STATES  ■  'fI3/ 

1  '  NO.  OF  PLANT  INPUTS  «  'f!3/ 

2  '  NO.  OF  DISTURBANCE  INPUTS  ■  'f  13/ 

4  '  NO.  OF  PLANT  OUTPUTS  »  'f  13/ 

5  '  NO.  OF  CONTROLLER  STATES  <EACH  CONTROLLER)  =  'fI3) 

WRITE(6f910) 

910  FORMAT < ;OPLANT  STATE  MATRIX  —  AP'> 

110  DO  112  I  *  IfNP 

READ(5f914)  <AF<IfJ)fJ=1fNP) 

112  WRITE<6f913)  < AP< I f J) f J=1 fNP) 

913  FORMAT  < '  'f8G13.6> 

914  FORMAT <6F 12.7) 

915  FORMAT  <  8G13 ♦ 6 ) 

WRITE <6 f 916) 

916  FORMAT ( "OPLANT  CONTROL  INPUT  MATRIX  —  B1P'> 

120  DO  122  I  ~  IfNP 

READ(5f914) (BIP(IfJ)fJ-IfNUP) 

122  WRITE<6f913X81P<IfJ>  fJ=IfNUP> 

WR!TE<6f918) 

918  FORMAT < 'OOBSERVER  MATRIX—  GE"  ) 

130  DO  132  I^IfNP 

READ<5f914) (GE( I f J) f J=1 fNWP) 

132  WRITE<6f913XGE<IfJ)fJ*1fNWP> 
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n  n  n  n  n 


CP'  ) 


WRITE  ( 6  x 920 ) 

920  FORMAT ( 'OPLANT  OUTPUT  MATRIX  -- 
140  HO  142  1=1 rNOP 

READ (5x  914 )  (CP( I x  J)  xJ=lxNF*> 

142  WRITE(6x913) (CP(IxJ) »J=1.NP) 

WRITE (6x 922) 

922  FORMAT ( 'OCONTROLLER  STATE  MATRIX  —  FC') 

ISO  DO  152  I  =  1  xNC 

READ ( 5  x914)(FC(IxJ) x J=lxNC) 

152  WRITE (6x913) (FC( I x J) x J=1 x  NC ) 

WRITE(6x924) 

924  FORMAT ( 'OCONTROLLER  CONTROL  INPUT  MATRIX  —  GC  ') 

160  DO  162  1=1 xNC 

READ( 5x914) (GC(Ix J) x J=lxNOP) 

162  WRITE(6x913)(GC(Ix  J)  xJ=lxNOF> 

WRITE (6x925) 

925  FORMAT ( 'OCONTROLLER  OUTPUT  MATRIX  (STATES)  —  HC') 

170  DO  172  1=1 xNUP 

READ (5x914)( HC (IxJ)xJ=lx  NC ) 

172  WRITE(6x913)(HC(IxJ) xJ=lxNC) 

WRITE(6x926) 

926  FORMAT ( 'OCONTROLLER  OUTPUT  MATRIX  (INPUTS)  —  EC') 

180  DO  182  1=1 xNUP 

READ (5x914) (EC(Ix J) xJ=lxNOP) 

182  WRITE ( 6 x 913 )  (EC  ( I  x  J )  x  J=1  xNQF') 

READ (5x928)  TxNT 
928  FORMAT (FI 0.4x15) 

XNT  =  NT 

DELTA  =  T/(XNT-1 ) 

WRITE(6x930)  TxNT 

930  FORMAT ( ' IT  *  'xF10.4/ 

1  '  NT  =  ' x 15/ 

3  '  T  =  SAMPLE  RATE.'/ 

6  '  DELTA  =  T/(NT-1)  =  INCREMENT  USED  IN  THE  NUMERICAL'/ 

7  '  INTEGRATIONS  TO  COMPUTE  PSITAUxPSITxPSIT2' / 

9  '  PSITAU1  USING  TRAPEZOIDAL  RULE.'//) 

WRITE (6x931) 

931  FORMAT (2Xx '  W  IS  THE  EXTERNAL  INPUT  (WHITE  GAUSSIAN  NOISE  WITH 

1  '  MEAN  =  O.Ox  AND  VARIANCE  *  1.0') 

GENERATE*WHITE*GAUSSIAN*NOISE*WITH*MEAN  =*0  AND%ARIANCE*=*i*** 

IX  -  11111 
DO  192  I  =  1x1200 
A  =  0.0 

DO  193  J  =  1x12 
CALL  RANDU< IXxIYxY) 

IX  =  IY 
193  A  -  A+Y 
192  W(I>  =  A-6 

DO  1199  I  =  IxNUP 
DO  1199  J  =  lxNP 
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ECCP(I»J>  =  0.0 
DO  1199  K  =  .1 1  NOP 

1199  ECCP(I.J)  =  ECCP(ItJ)  +  EC(I»K)*CP(KfJ) 

802  FORMAT (5X >5013. 6) 

835  FORMAT  <  5X  *  5013 « 6 ) 

C 

{ j  Jjt  jjt  ]jt  3)C  l([  )jt  )(t  $  )|(  )(l  DC  3^1  1C  )jc  ]f[  )j(  *  *  )j(  )j(  j(f  3^E  #  ijt  lj(  )j(  )^J  I(E  ]j( )([  1  3f(  lj(  )|[  )j(  lj(  )j(  )j(  *  )|i  ]j(  3|[  3|(  3jC  3f(  ^  i|t  )^E  )|t  ]|(  ])[  ](t  )|t !)( Ijt  ^  y  )|(  3(E  3^E  3|t  3)( 

C  CALCULATE  PH IT ( 0 >  > PHIT < DELTA )  > PHIT ( 2*DELTA ) > .  .  . >PHIT(T) 

tj  Ej(  l((  3)(  lj(  ijt  jjt  ]j(  #  3)C  3jC  *  3jC  3(c  *  )^E  *  I|E  3|C  $  3^  ijt  #  )(t  )fC  )|[  )j{  ]f!  I|(  )|[  J^C  3ft  3^!  3j^  3)C  ]|C  Jjt  3(C  3j(  30t  3^t  J^C  3jt  3f£  ]|C  )j(  3^C  3j(  3)C  3(C  3J(  3^C ))( 3|f  lf(  3(C  3jC  3jC  )J(  3jC  J^(  )j(  ]j[  3^C 

c 

DEL..HLF  ~  DELTA/2,0 
TI  ~  0,0 

DO  402  II  ==  1  >  NP 
DO  402  JJ  ~  1 >NP 
lFdl.EQ. JJ)  GO  TO  403 
PHIT1(II>JJ>1)  ~  0,0 
GO  TO  402 

403  PHIT1 ( 1 1  > JJ> 1 )  ~  1,0 

402  CONTINUE 

DO  4  II  -  2  >  NT 
TI  »:  TH  DELTA 
PHIT1 (lrlrll)  =  1, 

PHIT1 ( 1 >  2> 1 1 >  ==  (1,/10.)*(1.-EXP(-10.*TI)> 

PHIT1 (2>1>11)  =  0.0 
4  PHIT1 ( 2  >  2> 1 1 )  ~  EXP(-10.#TI ) 

DO  400  II  =  1 t NP 
DO  400  JJ  *  1 > NP 

400  PHXTdlrJJ)  "  PHIT 1  ( 1 1  > JJ > NT ) 

WRITE (6  > 860 ) 

860  FORMAT ( 5X  >  ' PHIT ' ) 

DO  361  I  =  1 > NP 

86.1  WRITE  ( 6  >  802 )  (PHIT ( I  >  J>  >  J~1  »NF*> 

N1  =•  0 

DO  1800  KK2  ~  1  >  6 
IAU  *  N1*0. 0125/5.0 
N1  ~  N1  +  1 
NTAU  ~  ( KK2-1 >#10 
NTAU  =  NTAU  +  1 

C  INITIAL  VALUE  YC2 ( -TIME+TAU ) >  YC3(-TIME)>  YC1 (-TIME) 

DO  31  I  »  1 > NUP 
YC2 ( I )  a  0.0 

31  YCKI)  ~  0.0 

DO  32  I  ~  1 > NUP 

32  Eld)  -  YCKI)  -  YC2(  I ) 

C  FROM  INITIAL  VALUE  XP(TIME)>  AND  YP(TIME)  ARE  EQUAL  TO  ZERO 

DO  35  1  =  1 >NP 

35  XP(I)  ■  0.0 

DO  36  I  ~  1 > NOP 

36  YP( I )  =  0.0 

C  FROM  INITIAL  VALUE  XCKTIME+T),  YC1 (TIME)  ARE  EQUAL  TO  ZERO 
DO  37  I  -  1  >NC 
XCl(I)  «  0,0 

37  XC2( I )  ■  0.0 

DO  52  I  «  1 > NUP 
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52  YCKI)  =  0.0 

C  ###########***#####***#***##**###*###***##*#***#*#*#*#*#******** 

c 

C  CALCULATE  PHIT(TAU) 

C 

I DEL  =  0 

DO  5000  I  =  If  NT 
I DEL  -  I DEL  +  1 
IF ( IDEL.EO.NTAU)  GO  TO  16 
GO  TO  5000 

16  DO  17  II  ~  1  r  NP 
DO  17  JJ  =  IfNP 

17  PHTAU  ( 1 1 f  JJ )  =  PHITKIIf  JJr  I  DEL) 

5000  CONTINUE 

C 

^  ^  ^  ^  ^  ^  j|^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

C  CALCULATE  PSIT(TAU) fFSIT<T> 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ** 
^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

C 

DO  550  I  ~  1 r NP 
DO  550  J  =  1 f NWP 

550  PS<IfJ)  *  0,0 
DO  551  1  =  1 » NP 
DO  551  J  a  If  NWP 
PSTAU(IfJ)  «  0.0 
PS IT  <  I  f  J)  =0.0 
DO  551  K  a  1 »NP 

PSTAU(IfJ)  a  PSTAU(IfJ)  +  PHIT1 < I fKfNTAU)#B1P<Kf  J ) 

551  PSIT(IfJ)  a  PSIT(IfJ)  +  PHIT1UfNfNT>*B1P(KfJ) 

DO  552  I  a  1 r NP 

DO  552  J  =  If NWP 

PSTAU(IfJ)  »  DELHLF*PSTAU(IfJ) 

552  PSIT(IfJ)  =  DELHLF#PSIT  < If J) 

60  DO  61  II  =  2fNT 
12  =  NT-11+1 

DO  62  I  =  1 f  NP 
DO  62  J  =  1 f  NWP 

62  PSIT(IfJ)  =  PSIT(IfJ)  +  PS(IfJ> 

DO  63  I  a  1 ,NP 

DO  63  J  =  l f  NWP 
PS(IfJ)  =  o.O 
DO  63  K  =  If NP 

63  PS(IfJ)  =  PS  < I f  J )  +  PHIT1<IfKfI2)*B1P(KfJ) 

DO  64  I  =  IfNP 

DO  64  J  =  1 f  NWP 

64  PS(IfJ)  =  DELHLF*PS<IfJ) 

DO  66  I  =  IfNP 

DO  66  J  a  1 fNWP 

66  PSIT<IfJ)  *  PS< I f J)  +  PSIT(IfJ) 

61  CONTINUE 

67  DO  68  I  =  IfNP 
DO  68  J  =  IfNWP 

68  PS  < I f  J )  =  0,0 
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IKINTAU.EQ.l)  GO  TO  55 
60  TO  69 

55  DO  58  IX  =  ItNP 
DO  58  JJ  “  1 r NWP 
58  PSTAU< II» JJ)  =  0.0 

GO  TO  77 

69  DO  70  II  =  2 » NTAU 
12  *  NTAU- I 1+1 

DO  71  I  =  1 r NP 
DO  71  J  =  1 » NWP 

71  PSTAU(IrJ)  =  PSTAU(ItJ)  +  PS(ItJ) 

DO  72  I  *  1 t NP 

DO  72  J  «  1 r NWP 
PS <  I » J )  -■  0,0 
DO  72  K  -  1 t NF 

72  PS<I»J)  “  PS(ItJ)  +  PHIT1(ItKtI2>#B1P(KtJ> 

DO  73  I  =  1 t NP 

DO  73  J  -  1 >  NWP 

73  PS(ItJ)  -  DELHl.F#PS ( I » J ) 

DO  76  I  ==  1  r  NP 

DO  76  J  =  1 t NWP 

76  PSTAU(ItJ)  -  F‘S<ItJ)  +  PSTAU<ItJ> 

70  CONTINUE 

77  CONTINUE 

C 

( ,  %  3(3  )4C  *  3(3  l|C  3(3  %  3(3  ]|k  *  3)3  3(3  l|>  3(3  3(3  lj(  3)3  3(3  ](t  3(3  l|(  3(3  3(3  ijt  l|(  3(3  I(t  3(1 3(3  3(3  3(3  3(3  3(3  3(3  (( 3(3  3(3  3(3  3(3  3(3  3(3  3(3 )(( 3(3  3)(  3(3  3(3  3(3  3(3  3(3  3(3  3(f  3(3  3(3  3(3  3(3  3(3  3(3  3(3  (3  3(3  Ij(  3(( 

C  START  TIME  LOOP 

TIME  =■•  0 

(j  (t  3(3  %  %  3(3  %  (( 3(3  *  *  %  )(3  3)3  lj(  3(3  %  %  ]( 3(3  %  3([  3|3  #  3(3  3(3  3(3  3(3  3(t  3(3  3(3  3(3  3(3  )(£  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3)3  3(3  3(3  3)3  3(3  3^  J(C 

c 

NN1  *  1 

DO  412  I  -  It  NOP 
412  UP(I)  ~  W(NN1 )~YCi < I ) 

l_  j  3(3  3(3  3(3  3f(  3((  3(3  3(3  j(3  3(3  3(3  3(3  )(3  3(3  3(3  3(3  ](3  !(3  3(3  3(3  ](3  )(3  )(3  3(3  j((  )(3  j((  j)3  3(3  3(3  j((  3)3  3(3  j(f  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3|C  3)3  3(3  3((  3(3  3(3  3(t  3(3  3(3  3(3  3(3  3(3  3(3  3(3  3fc  3(3  )(3  3)t  3(3  3(3  3(3  3(3  3(3  j(t  3(3 

C,  CALCULATE  XP ( TIME+TAU > t  AND  YP(TIME+TAU) 

C  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

DO  4"0  I  «  ItNUP 
430  E2 ( I )  -  YC1 < I )  -  YC2< I ) 

DO  220  III  »  It 400 

DO  253  1  »  ItNP 

0  »  0.0 

DO  254  J  -  ItNP 

254  Q  -  0  +  PHTAU<It J)#XP< J) 

253  XW1 < I )  b  0 

DO  255  I  b  lfNP 

Q  ®  0.0 

DO  256  J  b  1 tNUP 

256  0  =  0  +  PST  AU ( I t  J  >  *UP ( J  > 

255  XW2< I )  b  Q 

DO  257  I  s  ItNP 

257  XPTAU(I)  b  XW1 < I )  +  XW2(I) 

DO  260  I  b  1 tNOP 

G  b  0,0 

DO  261  J  «  ItNP 
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261  a  =  0  +  CPd»J>*XPTAU< J) 

260  YPTAU(I)  =  Q 

C,  He  He  )j(  )|( )(( He  He  )jc  He  He  i([  )jt  ]j(  )jt  He  i)(  ijt  )|(  ]jc  )j(  i(t  He  )|{  )|f  !j(  He  He  ijt  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  H 

C  TIME  =  TIME  +  T 

C  CALCULATE  XC2 ( TIME  +  T+TAIJ )  *  AND  YC2<  TIME+TAU) 

DO  300  I  ~  1 t NUP 
Q  ~  0,0 

DO  301  J  -  1  »NC 

301  Q  -  Q  +  HCd»J>*XC2(  J) 

300  XUIKI)  =  0 

DO  302  I  ==  1  r  NUP 

0  »  0.0 

DO  303  J  ~  1 t NOP 

303  Q  ~  Q  +  EC<Ir J)HeYPTAU(J) 

302  XW2 ( I )  ~  0 

DO  30-1  1  «  1 1  NUP 

304  YC2  < I )  --  XWl(I)  +  XW2(  I ) 

DO  305  I  ~  1 t NC 

Q  •-  0.0 

DO  306  J  »  1»NC 

306  0  -^0  +  FC  <  I » J  )  HeXC2  ( J  ) 

305  XWld)  ■  Q 

DO  307  I  ~  1 r  NC 

Q  »  0.0 

DO  308  J  »  1 » NOP 

308  0  ^0  +  GC(If  J)HeYPTAU(J) 

307  XW2  < I )  *  Q 

DO  309  I  -  1 »NC 

309  XC2  ( I )  ®  XWld)  +  XW2CI) 

L  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  H  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  H  He  He  He  He  He  He  He  H^  He  He 

C  CALCULATE  El 

DO  290  I  =  1 i NUP 
290  Eld)  “  YCl(I)  -  YC2(  I ) 

(.i  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He 

C  CALCULATE  XP(TIME)>  AND  YP(TIME) 

1/  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  H  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  H 

373  DO  500  I  ~  1  »NP 

Q  =  0.0 

DO  501  J  ~  1 t NP 

501  a  ~  a  +  PHITdr J)HtXP(J) 

500  xwid)  =  a 

DO  502  I  -  1 f NP 

a  -  o.o 

DO  503  J  “  1 t NUP 

503  a  «  a  +  PsiTdr  j>*up<  j) 

502  XW2( I )  »  Q 

DO  504  I  =  1 f NP 

504  XP(I)  =  XWl(I)  +  XW2  < I ) 

DO  507  I  *  1 t NOP 

a  =  o.o 

DO  508  J  -  1 » NP 
508  Q  "  Q  +  CPdf J)*XP( J) 
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507  YP(I)  -  Q 

C  #  ###*######*####***  #**##***##*##  **%***%*****  *********  ******-t**.*#  + 

C  CALCULATE  XCK2*TIME)f  AND  YC1( TIME  > 

C  itc*####*******#*###*#*  i**######)^####***#########********  ****** 

no  700  I  ■=  1 1 NUP 

Q  •=  0.0 

no  701  J  ”  1 f NC 

701  a  =  Q  +  HC<I»J>*XC1(J> 

700  XWKI)  ~  Q 

no  702  I  ~~  If  NUP 

n  =  o.o 

no  703  J  =  1 f NOP 

703  Q  =  0  +  EC( I »  J>*YP< J) 

702  X W2 ( I )  ~  0 

no  704  I  -  If NUP 

704  YCl(I)  =  XWKI)  +  XW2(  I ) 
no  705  I  ~  IfNC 

0  ~  0.0 

no  706  J  -  IfNC 

706  a  -  a  +  fc<IfJ)*xckj> 

705  XWKI)  =  Q 

no  707  I  ~  IfNC 

Q  ~  0.0 

DO  708  J  "•  I  f  NOP 

708  U  ~  Q  +  GC<  I  f  ,J)*YP(  ,J) 

707  XW2 ( I )  ~  0 

no  709  I  =-  IfNC 

709  XCKI)  =  XWKI)  +  XW2  <  I ) 

NN1  “  NN1  +  1 

DO  221  I  =■  IfNOP 
221  UP(I)  ~  W(NNl)-YCKI) 

tj  ■(t  ^  t  3(1  ^  3(1  ^  3(1  ^  ^  3|[  )j(  3(1  l(f  i(C  l|c  3(1  3(1  #  )((  ]|(  ](( 3(1 ))( )(t  )|t  3(1  )j{  3(1  3(1  $  $  3(1  )j(  )|!  )(t  ](t  3(1  )jl  ]|C  3(1  3(1  3^!  ])!  I(t  3^  3(1  ]|C  3(1  3(1  3(1 )(( )(( 3(1 

C  CALCULATE  E2 

(„■  3(1  #  *  jjc  l(t  3(1  l(t  3(3  3(1  )j{  if!  3((  %  3(1  3(1  )|(  ])[  ijt  3#C  l(!  3(1  ]j(  3(1  3j(  3(1  Ij(  ]f(  ]j(  1(1 3(1 J(l  3(1 3(1 3(1 3(1 3(1 3(1 3((  3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1  )(t  3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1 3(1 3^  3(1 3(1 3(1 3(l 

no  540  I  =  If NUP 
540  E2(I)  »  YCl(I)  -  YC2(I) 

HO  560  I  ~  If NUP 
EA(IIl)  *  EKI) 

560  EB(IIl)  “  E2<I) 

220  CONTINUE 

l 

L  CALCULATE  THE  STEADY  STATE  SAMPLE  COVARIANCE  OF  ERRORS 

( 

SMEANA  0.0 
SMEANB  =--0.0 
1*0  561  I  501  f  800 
.MEANA  *  SMEANA  +  EA(I) 
rtlANB  SMEANB  +  EB  ( I ) 

ME ANA  SMEANA/300 
s  m  AN  it  SMEANB/300 

A  A  0.0 
«M*  0.0 
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HO  563  I  =  501 f 800 
VEA  -  VEA  +  < EA ( I ) -SMEANA ) *# 

563  UEB  *  UEB  +  <EB< I )-SMEANB>** 
PEASS(KK2)  -  VEA/300 
F‘EBSS(KK2 )  -  UEB/300 

1800  CONTINUE 
WRITE! 6 » 570)  TAU 

570  FORMAT <  5X » ' PEASS (  TAU  =  '»F12.8»') 
WRITE (6 f 571) (PEASS(I) *1-1,6) 

571  FORMAT <5XrF 18. 10) 

WRITE ( 6 r 572 )  TAU 

572  FORMAT < 5X » ' PEBSS<  TAU  *  '»F12.8»') 
WRITE ( 6  >  571 ) <  F'EBSS  (I)»l  =  l?6) 

1801  CONTINUE 
STOP 
END 

SUBROUTINE  RANDU < I X  *  I Y i YFL ) 

IY  -  IX*65539 
IF<IY)5»6,6 

5  IY  *  IY  +  2147483647+1 

6  YFL  «  IY 

YFL  ■"  YFL#0 . 46566 13E-9 

RETURN 

END 
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ro  rj 


EXAMPLE - 27 H  ORDER  PLANT,  1ST  ORDER  CONTROLLER 

NO.  OF  PLANT  STATES  ~  2 
NO.  OF  PLANT  INPUTS  =  1 
NO.  OF  EXTERNAL  INPUT  =  J 
NO.  OF  PLANT  OUTPUTS  =  1 

NO.  OF  CONTROLLER  STATES  (  EACH  CONTROLLER)  =  1 


PLANT  STATE  MATRIX  —  AP 

♦  0  1.0 

♦0  -10.0 

PLANT  CONTROL  INPUT  MATRIX  —  Bp 

.0 

200.0 


PLANT  OUTPUT  MATRIX  —  CP 

1.0  .0 


CONTROLLER  STATE  MATRIX  —  FC 
.523810 

CONTROLLER  CONTROL  INPUT  MATRIX  —  GC 
- • 18162 

CONTROLLER  OUTPUT  MATRIX  (STATES)  —  HC 

1.0 


CONTROLLER  OUTPUT  MATRIX  (INPUTS)  —  EC 
1.381 
NT  =■  51 

7  =  SAMPLE  PERIOD  =  0.0125  SEC 

DELTA  *  T /  (NT-1 )  »  INCREMENT  USED  IN  THE.  NUMERICAL 

INTEGRATIONS  TO  COMPUTE  PSITAU»PBIT»PSIT2» 
PSITaUI  USING  TRAPEZOIDAL  RULE. 

W  IS  THE  EXTERNAL  INPUT  (WHITE  GAUSSIAN  NOISE  WITH 
MEAN  -  0.0,  AND  VARIANCE  -  1.0 

7 HE  STEADY  STATE  SAMPLE  VARIANCE  OF  ERRORS 
PEASS 

0.0 

0.0002195683 

0,0008547062 

0,0018880414 

0.0033245913 

0.0051899776 

PEBSS 

0.0051899776 

0,0033477221 

0,0019137899 

0.0008709673 

0.0002244494 

0.0 
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; '  «•»■' l."1 


APPENDIX  C 


DESIGN  OF  THE  STATE  OBSERVER 

In  section  III,  y  (tw)  the  Input  to  channel  1,  go  through  the 
blocks  named  1st  DIGITAL  CONTROLLER  AND  OBSERVER.  The  output  of 
the  first  block  serves  as  the  input  to  the  plant  and  the  output  of  the 
second  is  the  estimate  of  the  state  variables  of  the  plant  and  is  used 
to  estimate  y£2(tk  +  **).  This  appendix  describes  how  to  estimate  the 
state  variables  XpU^)  from  yp( tk) - 

1 .  DESIGN  OF  STATE  OBSERVER 

Reference  3  defines  an  observer  as  the  subsystem  that  estimates 
the  state  variables  of  a  dynamical  system,  based  on  measurements  of 
the  input  up(t)  and  the  output  yp(t).  Figure  22  shows  the  block 
diagram  of  an  observer  which  is  formulated  as  a  feedback  control  with 
Gc  as  the  feedback  matrix.  The  design  objective  is  to  select  the 
feedback  matrix  Ge  such  that  yp(t),  the  estimate  of  yp(t),  will  approach 
yp(t)  as  fast  as  possible.  When  yp(t)  equals  yp(t),  the  dynamics 
of  the  state  observer  are  described  by 

x*(t)  =  ApXp(t)  +  Bpup(t)  C-l 

which  Is  identical  to  the  state  equation  of  the  system  (plant)  to 
be  observed.  In  general,  with  up(t)  and  yp(t)  as  inputs  to  the 
observer,  the  dynamics  of  the  observer  are  represented  by 

xj(t)  «  [Ap  -  GeCp]xp(t)  +  BpUp(t)  +  Geyp(t)  C-2 
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Since  yp(t)  equals  Cpxp(t),  the  equation  C-2  Is  written  as 

x*(t)  =  Apxp(t)  +  Bpup(t)  +  GeCp[xp(t)  -  xf(t)]  C-3 

The  significance  of  this  expression  is  that  if  the  initial  values 
of  xp(t)  and  xp(t)  are  identical,  the  equation  reverts  to  that  of 
equation  C-l,  and  the  response  of  the  observer  will  be  identical  to 
that  of  the  original  system.  [In  the  model  in  Section  III,  the  initial 
value  of  x*(tk)  is  unknown.  Therefore,  the  design  of  the  feedback 
matrix  Gfi  for  the  observer  is  significant  only  if  the  initial  conditions 
of  xp(t)  and  xp(t)  are  different. 

If  we  subtract  equations  (C-3)  from  (C-l),  we  have 

[xp(t)  -  x*(t)]  «  [Ap  -  GeCp][xp(t)  -  xj(t)]  C-4 

which  may  be  regarded  as  the  homogeneous  state  equation  of  a  linear 
system  with  the  coefficient  matrix  [A  -  GeCp] .  The  characteristic 
equation  of  [A  -  GfiCp]  and  of  the  state  observer  is  then 

| X  -  (A  -  GeCp) |  =  0  C-5 

Since  we  are  interested  In  driving  xp(t)  as  close  to  xp(t)  as  possible, 
the  objective  of  the  observer  design  may  be  stated  as  to  select  the 
elements  of  Gg  so  that  the  natural  response  of  equation  C-4  decays  to 
zero  as  quickly  as  possible.  In  other  words,  the  eigenvalues  of 
[A  -  GeCp]  should  be  selected  so  that  xp(t)  approaches  xp(t)  rapidly. 
However,  it  must  be  kept  In  mind  that  the  approach  of  assigning  the 
eigenvalues  of  [A  -  GeCp]  may  not  always  be  satisfactory  for  the 
purpose  of  matching  all  the  observed  states  to  the  real  state,  since 
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the  eigenvalues  control  only  the  denominator  or  polynomial  of  the 
transfer  relation,  while  the  numerator  polynomial  is  not  controlled. 
More  details  and  an  example  of  the  statement  above  are  available  in 
Reference  4. 

2.  EXAMPLE 

The  following  example,  which  is  the  example  used  in  Section  II 
and  III,  is  used  to  illustrate  the  technique  described  above. 

From  the  example  in  Section  II,  we  have 

Ap  =  0  1  C-6 

0  -10 

Bp  =  0  C-7 

200 

and 

Cp  *|l  °|  C-8 

Let  the  feedback  matrix  be  designated  as 

Ge  ■  9e)  C-9 

9,2 

Substitution  of  equations  (C-6),  (C-8),  and  (C-9)  into  (C-5),  gives 


Then  the  characteristic  equation  of  the  state  observer  is 


Xx2  +  (10  +  ge1)x  +  (10gel  +  ge2)  -  0  C-11 

Let  the  eigenvalues  of  XI  -  (Ap  -  GeCp)  be  x  =  -15,  -15  then 
the  characteristic  equation  should  be 

X2  +  30x  +  225  =  0  C-12 


Equating  like  terms  of  equations  (C-11)  and  (C-12)  gives 


9el  =  20 
ge2  *  25 


Figure  23  illustrates  the  responses  xpl(t)  and  xpl(t)  for  the 
following  initial  states 


xp(0>  * 


0  , 
0 


2 

0 


Shown  in  the  same  figure  is  the  response  of  xpj(t)  when  the  state 
observer  is  designed  for  eigenvalues  at  x  =  -20,  -20;  in  this  case 
gel  *  30  and  ge2  *  100.  However,  it  is  seen  from  the  figure  that 
xpl  for  both  gel  =  20,  ge2  =  25  and  gel  =  30,  ge2  *  100  are  approxi¬ 
mately  the  same  deviation  from  xp. 

Figure  24  Illustrates  the  response  xp2(t)  and  xp2(t)  for  the 
two  cases  of  observer  design.  The  characteristics  of  xp2  for  both 
gel  »  SO,  ge2  ■  25  and  ge1  *  30,  ge2  *  100  are  the  same  as  that  of 
both  x*j  which  are  approximately  the  same  deviation  from  xp2. 

As  mentioned  earlier,  the  eigenvalues  may  not  always  be  satis- 
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FIGURE  23  :  STATE  xlp  AND  THE  OBSERVED  STATES 
OF  THE  OBSERVER 


factory  for  the  purpose  of  matching  all  the  observed  states  to  the 
real  states.  Reference  4  shows  that  the  selecting  larger  values  for 
gej  and  ge£  to  give  faster  transient  response  for  the  observer  is 
not  always  best.  Sometimes,  the  large  values  of  gej  and  ge2  will 
only  give  faster  transient  response  for  one  of  the  observed  states 
but  not  the  other  and  the  details  are  in  Reference  4. 

Since  both  xpj  and  x^  converge  to  xp^  and  xp2  respectively  at 
the  same  rate  for  x  =  -20,  -20  and  x  =  -15,  -15,  then  we  can  select 
either  set  of  eigenvalues  for  this  example. 
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APPENDIX  D 

SOFTWARE  FOR  THE  NEW  MODEL 

1.  FLOWCHART  AND  DESCRIPTION  OF  MAJOR  COMPONENTS  AND  SUBROUTINE 
The  main  program  in  the  software  for  the  new  model  is  called 
PROGRAM  SKEW1 .  Its  major  computational  tasks  are  to  develop  the  state 
variable  model  of  the  complete  closed-loop  system  and  to  compute  the 
steady-state  covariance  of  the  states,  the  controller  outputs  and  the 
errors  eA  and  eg.  This  program  is  almost  the  same  as  PROGRAM  SKEW, 
described  in  Appendix  A,  except  that  the  inherent  errors  are  computed 
from  the  output  of  channel  2,  y^f^ic  +  t),  and  its  estimated  value 
from  the  input  of  channel  1,  y£2(tk  +  **)•  All  the  variables  from 
Appendix  A  plus  the  variables  for  calculating  the  estimated  value  of 
the  output  of  the  output  of  channel  2  are  used  in  this  proqram. 

A  flowchart  for  this  program  appears  in  Figure  25.  The  first 
block  shows  the  data  input.  The  variables  are  self-explanatory  and 
discussed  in  Appendix  A  except  for  a  new  variable  NTAU1 ,  which  represents 
the  value  of  t*  as 


x*  =  (NTAUli  -  D*  T 

trtt : i) 

and 

NTAUli  =  ITAU1  *  NT  =  1 

RTF 

for  ITAU1  *  0,  1,  2 . NT2 

where  NT  must  be  an  integer. 

NT? 


D-l 


D-2 


In  block  two,  gaussian  noise  is  computed  by  using  subroutine  RANDU. 
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Read  in  heading  and  description  of  the  plant  and 
the  controller ;np ,nup .n^ ,nQp ,nc , Ap . Blp , Cp , Ec , Fc , 

Gc|Hc,Ge>T|NT‘  811(1  m2 _ 


Compute  the  Gaussian  white  noise  by  using 
SUBROUTINE  RANDU 


Calculate  <Kt*)  ,  and  $(T)  by  using  eq.2-4 


Calculate  'Kt;L),  ,  YCT),  and  ^(T)  by  using 

the  trapezoidal  rule  for  numerical  integration 


Calculate  x*(tR) ,  x*(tR+T*) ,  y*(tR+T*) ,  x*2 (tR+1+t*) , 
and  yc2^tK+T*)  by  using  eqs.  3-13,  3-15,  3-16,  3-7, 
and  3-8 


Calculate  xc2 ^fcK+l+T^  ’  yc2^tK+T^  by  usin8  e9s.  3-9, 
and  3-10 


Calculate  xc^(tR+^),  yc^(tR)  by  using  eqs.  3-5,  and  3-6 


Calculate  eA,eB,PeAgs,  and  PeBgg  by  using  eqs.  3-17, 
3-18,  3-19,  and  3-20 


FIGURE  25  :  FLOWCHART  DESCRIBING  THE  MAJOR  COMPUTATIONS 
OF  PROGRAM  SKEW1 


The  third  block  specifies  the  calculation  of 
and  ${T). 

The  fourth  block  specifies  the  calculation  of  ^(t^, 
and  i|i(T)  by  using  the  same  technique  in  Appendix  A. 

In  block  five,  xj(tk).  xp(tk  +  x*),  x*2(tk+1  +  x*),  and 
yc2(tk  +  x*)  are  calculated  from  equations 

xp(tk)  ~  #(tk,  tk_] )xp(tk_i )  +  tHtk,  tk_ i ) Up ( tk_ i ) 

+  ^(tk.  tk-i  )yp(tk_i ) 

xp(tk  +  T*)  =  ♦(tk  +  X*.  tk_i )x*(tk_1 )  +  *(tk  +  X*.  tk_i )Up(tk_1 ) 
+  4>l(tk  +  X*.  tk_i )yp(tk_i ) 

xc2^k+l  +  T  )  =  Fcxc2^k  +  T*^  +  ®c^p(^k  +  T*) 
yc2(tk  +  T*)  3  Hcxc2(tk  +  T*)  +  Ecyp(tk  +  T*) 

In  block  six,  xc2(tk+j  +  x)  and  yc2(tk  +  t)  are  calculated 
from  equations 

xc2^k+l  +  =  Fcxc2^k  +  T)  +  typttfc  +  T) 

yc2(tk  +  t)  =  Hcxc2(tk  +  x)  +  Ecyp(tk  +  x) 

In  block  seven,  xcj(tk+])  and  yC](tk)  are  calculated  from 
equations 

xcl(W  *  Fc»cl  <*k>  +  Vpftk) 
yc,(tk)  ■  Hcxc1(tk)  +  Ecyp(tk) 
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MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS- 1963-A 


In  block  eight,  the  inherent  errors,  e^,  eB,  and  the  steady- 
state  covariance  of  errors:  PEASS,  PEBSS,  are  computed  from  equations 


eA^) =  yci<V  ■ 

tk  +  r  <  t  <  tk+1,  0  <  t  <  T,  k  =  0,  1,  .  .  . 
eB(t)  =  -  yC2^k +  t) 


tk+i  K  tfc+i  +  T,  0  <  T  <  T,  k  =  0,  1 ,  . 


N 


PEASS  =  1 


N  i=201 
N 

PEBSS  =  1 

N  i=201 


^eAi  "  ffA^  ^eAi 
teBi  "  ®B^  ^eBi 


-  *a3T 


V 


Since  the  majority  of  this  program  is  the  same  as  the  program  in 
Appendix  A,  then  only  new  arrays  will  be  shown  in  the  following 
XP1  (NP)  XC3  (NC) 

YP1  (NOP)  XPTAU1  (UP) 

YC3  (NUP)  PSIT2  (NP,  NUP) 


2.  EXAMPLE 

The  computer  listing  for  PROGRAM  SKEW1,  of  the  new  model,  written 
in  FORTRAN,  appears  in  Appendix  E.  For  each  value  of  ,  PEASS  and 
PEBSS  for  each  value  of  t*  are  shown  at  the  end  of  the  listing.  To 
compute  PEASS  and  PEBSS  of  the  new  model,  the  system  has  been  waiting 
until  the  state  observer  xp  is  equal  or  close  to  xp  and  the  system  is 
in  steady-state  transition. 

From  Appendix  C,  the  state  observer  xp  equals  xp  at  approximately 
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30T  (gej  *  20,  ge2  *  25).  Thus,  PEASS  and  PEBSS  can  be  computed  at 

the  same  time  (between  200T  and  300T)  as  PEASS  and  PEBSS  of  the  basic  { 

.  i 

model  In  Appendix  A  are  computed. 

i 

( 

i 

1 

1 

1 

I 

\ 

I 

I 

I 

I 

I 

I 

n 

fi 
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APPENDIX  E 

COMPUTER  PROGRAM  LISTING 
FOR 

PROGRAM  SKEW!  AND  EXAMPLE 
OF  OUTPUT  WRITTEN  IN  FORTRAN 
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CCCCC  PROVIDE  MAXIMA  FOR  CALLED  ARRAYS 
NPM  =  2 

NUPM  =  2 
NWPM  «  2 
NORM  *  1 

NCM  *  2 

NHM  =  NPM  +  NUPM 
NFM  “  NPM  +  2*NCM 
NRRM  ~  2*NPM  +4 

C 

C  READ  INPUT  DATA 


WRITE (6t 899 ) 

89V  FORMAT ('I') 

100  READ<5»900)  ID 
900  FORMAT (20A4) 

WRITE <6» 902)  ID 
902  FORMAT (  '  1 ' »20A4 ) 

READ <  5 » 906 ) NP , NUP » NWP » NOP  r  NC 
906  FORMAT (51 3) 

WRITE <6, 908 )  NP » NUP , NUP » NOP » NC 
908  FORMAT  < 'ONO,  OF  PLANT  STATES  ■  '*13/ 

1  '  NO.  OF  PLANT  INPUTS  =  ,13/ 

2  '  NO,  OF  DISTURBANCE  INPUTS  *  '»  13/ 

4  '  NO,  OF  PLANT  OUTPUTS  ■  13/ 

5  '  NO,  OF  CONTROLLER  STATES  (EACH  CONTROLLER)  =  ',13) 

WRITE(6»910) 

910  FORMAT ( 'OPLANT  STATE  MATRIX  —  AP') 

110  DO  112  I  =  1 »NP 

READ ( 5 f 91 4)  ( AP( I » J) » J=1 » NP) 

112  WRITE(6»913)  ( AP( I » J ) » J=1 , NP ) 

913  FORMAT ( '  '»8G13,6) 

914  FORMAT (6F12,7) 

915  FORMAT (8G13. 6) 

WRITE(6»916) 

916  FORMAT < ' OPLANT  CONTROL  INPUT  MATRIX  —  B1P') 

120  DO  122  I  -  1 »NP 

READ(5»  914)  (BlP(IrJ)i>J=l  »NUP) 

122  WRITE(6»913XB1P(I»J) > J=1*NUP) 

WRITE(6»918) 

918  FORMAT (  ' OOBSERVER  MATRIX—  GE'> 

130  DO  132  1=1 >NP 

READ (5, 91 4) (GE( I r J) r  J=1 »NWP) 

132  WRITE(6»913XGE(IfJ)fJ=l  »NWP) 

WRITE(6f 920) 

920  FORMAT < 'OPLANT  OUTPUT  MATRIX  —  CP') 

140  DO  142  1=1 »NOP 

READ(5»914) <CP(I»J)»J=1*NP) 

142  WRITE<6»913XCP(I»J)»J=1,NP> 

WRITE(6»922) 

922  FORMAT < 'OCONTROLLER  STATE  MATRIX  —  FC') 

150  DO  152  I  »1»NC 
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o  u  u 


152 


READ(5»914)(FC<I»J)f J=1»NC) 

WRITE<6f913HFC<I»J>f  J-lfNC) 

WRITE (At  924) 

924  FORMAT < 'OCONTROLLER  CONTROL  INPUT  MATRIX  —  GC  ') 

160  HO  162  1 1  f  NC 

READ(5f914HGC(If J)f J=lfNOP> 

162  WRITE (6  f  913) ( GC ( I f  J ) f  J= 1 1 NOP ) 

WRITE (6 .925) 

925  FORMAT ( 'OCONTROLLER  OUTPUT  MATRIX  (STATES)  —  HC') 

170  DO  172  l~lt NUP 

READ <5 *91 4)  (HC( I f J) f J~1 f NC) 

172  WRITE (6.913) (HC(If J) r J=1*NC) 

WRITE(6f 926) 

926  FORMAT < 'OCONTROLLER  OUTPUT  MATRIX  (INPUTS)  —  EC') 

18 0  DO  182  I~1  f  NUF' 

READ(5f914> (EC<If J) f J=lfNOP) 

182  URITE(6»913) (EC ( I r J) r J=1 fNOP ) 

READ (5.928)  T f  NT 
928  FORMAT (F10.4fI5> 

XNT  *  NT 

DELTA  ~  T/(XNT-1 ) 

WRITE (6 » 930)  T.NT 

930  FORMAT ( ' IT  »  'r FI 0.4/ 

1  '  NT  »  '.15/ 

3  '  T  *  SAMPLE  RATE.'/ 

6  '  DELTA  »  T/(NT-1)  *  INCREMENT  USED  IN  THE  NUMERICAL'/ 

7  '  INTEGRATIONS  TO  COMPUTE  PSITAU. PSIT f PSIT2' / 

9  '  PSITAU1  USING  TRAPEZOIDAL  RULE.'//) 

WRITE ( 6.931) 

931  FORMAT < 2X » '  W  IS  THE  EXTERNAL  INPUT  (WHITE  GAUSSIAN  NOISE  WITH 

1  '  MEAN  *  0,0.  AND  VARIANCE  *  1.0') 

C 

C  #****#*#r**##JF*****##********#****#*###**##***#*#**#**# ********** 

C  GENERATE  WHITE  GAUSSIAN  NOISE  WITH  MEAN  »  0  AND  VARIANCE  »  1 
C  ***************************************************************** 

c 

ix  »  mu 

DO  192  I  ~  If  1200 
A  a  0.0 

DO  193  J  “  lr 12 
CALL  RANDUdXf  IYfY) 

IX  *  IY 
193  A  -  A+Y 
192  W(I)  *  A -6 

DO  1199  I  *  If NUP 
DO  1199  J  IfNP 
ECCP(IfJ)  «  0.0 
DO  1199  N  =  1 r  NOP 

1199  ECCP ( I *J>  -  ECCP(IfJ)  +  EC<IfK)*CP(K»J) 

802  F0RMAT(5Xf5Gl3.6) 

835  FORMAT (5Xf5G13.6) 

***************************************************************** 

CALCULATE  PHIT(0>»PHIT (DELTA) »PHIT(2*DELTA) ,  .  .  . fPHIT(T) 
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Ci  ^  ^  ^  ^|t  He  He  1 3|(  He  He  J  He  *  He  *  i|t  i|[  ]|(  ]|t  He  ijc  )j(  i|[  i|(  He  He  H  He  He  He  He  H  He  He  H  He  He  H  H  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He  He 
C 

DELHLF  =  DELTA/2.0 
TI  «  o.O 

DO  402  II  =  1 pNP 
DO  402  JJ  =  1 f  NP 
IF(II.EO.JJ)  GO  TO  403 
PHIT1  (IlfJJrl)  =  0.0 
GO  TO  402 

403  PH I T 1  ( 1 1  f  J J » 1 )  =  1.0 

402  CONTINUE 

DO  4  11  =  2» NT 
TI  =  TI+DELTA 
PHITKlf  lfll)  =  i, 

PHITl(lr2»Il)  *  (1./10. >*<1.-EXP<-10.*TI)  ) 

PHIT1 <2.1.11)  =  0.0 
4  PHIT1  (2.2.11)  =  EXP<-10.HeTI  > 

DO  400  II  --  1 .  NP 
DO  400  JJ  *  1 .  NP 

400  PHIT< I I » JJ)  ■  PHITKII.JJ.NT) 

N1  =  0 

DO  1801  KK1  ■  1.6 
NTAU1  *  (KKl-l)HelO 
NTAU1  -  NTAU1  +  1 
TAUU  »  NIHeO. 0125/5. 

N1  =  N1  +  1 
DO  1800  KK2  ■  1.6 
NTAU  =  (KK2-1 ) *10 
NTAU  -  NTAU  +  1 

C  INITIAL  VALUE  YC2< -TIME+TAU) .  YC3<-TIME).  YCl(-TIME) 

DO  31  I  =  1 .  NUP 
YC2( I )  ■  0.0 
YC3( I )  »  0.0 

31  YC1 < I )  =  0.0 

DO  32  I  =  l.NUP 

32  El ( I )  =  YC3< I )  -  YC2 < I ) 

C  FROM  INITIAL  VALUE  XP(TIME),  AND  YP(TIME)  ARE  EQUAL  TO  ZERO 
DO  35  I  =  1 .  NP 
XP1 < I )  a  2.0 

35  XP<I)  a  0.0 

DO  435  I  »  1? NOP 

a  =  o.o 

DO  436  J  *  1 »NP 
436  Q  a  q  +  cp(I» J>*XP1<J> 

435  YP1 < I >  a  a 

DO  36  I  =  lr NOP 

36  YP(I)  a  o.O 

C  FROM  INITIAL  VALUE  XCl<TIME+T).  YCl(TIME)  ARE  EQUAL  TO  ZERO 
DO  37  I  =  1 f  NC 

XC1<I>  a  0.0 
XC3< I >  a  0.0 

37  XC2( I )  »  o.O 

DO  52  I  =  IfNUP 
52  YCUIt  »  0.0 
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c 

c  calculate:  phit<tau> 

c 

L  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c 

I DEL  -  0 

DO  5000  1  *  lr NT 
I DEL  =  I DEL  +  1 

1 F  <  NT  AU . EQ . NTAU1 « AND . I DEL . EQ « NTAU )  GO  TO  7 
IF( IDEL » EQ .NTAU)  GO  TO  16 
T F ( I DEL . EQ ♦ NT AU 1 )  GO  TO  22 
GO  TO  5000 

7  DO  8  II  -  KNP 
DO  8  JJ  -  IfNP 

PHTAU(IIfJJ)  «  PHITKIIf JJ»IDEL) 

8  PHTAUKIIf  JJ)  ~  PHIT1  < II f JJf  I  DEL) 

GO  TO  5000 

16  DO  17  II  -  IfNP 
DO  17  JJ  ~  l.NP 

17  PHTAU(IIfJJ)  ~  PHZTKIIr  JJf  I  DEL) 

GO  TO  5000 

22  DO  23  II  »  IfNP 
DO  23  JJ  “  IfNP 

23  PHTAUKIIf  JJ)  “  PHITKIIf  JJf  IDEL) 

5000  CONTINUE 

C 

C  CALCULATE  PSIT(TAU) fPSIT(T) 

f  4k  k  k  Ik  da  ik  -Ja  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  aka  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k 

T  f  f  ™  k  T  T  T  T  T  f  f  ^  k  k  f  f  ^  f  k  ^  ^  k  k  k  k  k  k  k  ^  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k 

c 

DO  550  I  ~  ltNP 
DO  550  J  «  IfNWP 
PKIfJ)  =  0.0 

550  PS(IfJ)  «  0.0 

DO  551  I  =  IfNP 

DO  551  J  IfNWP 

PSIT2< I f J)  »  0.0  • 

PSTAU(IfJ)  -  0.0 
PSTAUl(IfJ)  =  0.0 
PSIT< I f J)  «  0.0 

DO  551  K  -  IfNP 

PSIT2(  I  f  ,J)  *  PSIT2(IfJ>  +  PHITl(IfKfNT)4cGE(Kf  J) 

PSTAU(IfJ)  *  PSTAU(IfJ)  +  PHIT1  <  I  r  Kf  NTAU)4cBlP(Kf  J) 

PSTAUKIfJ)  -  PSTAUKIfJ)  +  PHIT1  ( I  f  KrNTAUl  )4cBlP(Kf  J) 

551  PSIT(IfJ)  *  PSIT(IfJ)  +  PHITl<IfKrNT)4iBlP<Kf J) 

DO  552  I  «  IfNP 

DO  552  J  -  IfNWP 

PSIT2(  I  f  J)  *  DELHLF4(PSIT2<If  J) 

PSTAU(IfJ)  =  DELHLF4(PSTAU(If  J) 

PSTAUKIfJ)  =  DELHLF4tPSTAUl  <  I » J) 

552  PSIT  <  I  f  J)  -  DEl.HLF4iPSIT<If  J) 

60  DO  61  II  ~  2f NT 

12  *  NT-Il-f  1 
DO  62  I  ■  IfNP 
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62 


63 


64 


o>6 

61 


80 

85 


86 


90 


6"> 


93 


94 


96 

91 

67 

68 


58 


DO  62  J  =  1 t NWP 

F'SIT2<  I  f  J>  =  PSIT2<  I »  J)  +  P1(I,J) 

PSIT(IfJ)  =  PSIKItJ)  +  PS(IrJ) 

DO  63  I  -  1,NP 
DO  63  J  -  1 t NWP 

F-Klrj)  a  0.0 

PS(I,J)  a  0.0 
DO  63  K  =  1 »NP 

P1(X,J)  a  PKIrJ)  +  PHIT1  ( I»  K»  I2)*GE(Kr J) 
PS(I,J)  a  PS(IfJ)  +  PHIT1 ( I»K»I2)#B1P(K»J) 
DO  64  I  =  1 »NP 

DO  64  J  a  1 f NWP 
PKIrJ)  a  DELHLF*Pl(lf J) 

PS(IfJ)  a  DELHLF*PS<If J) 

DO  66  1  a  i,NP 
DO  66  J  «  1 t NWP 

PSIT2( I » J )  a  P1<I,J)  +  PSIT2 ( I >  J ) 

PSIT< I f J)  a  PS(I»J)  +  PSIT<I»J) 

CONTINUE 

DO  80  I  «  1  ,NP 

DO  80  J  a  1 T NWP 

PS(IrJ)  a  o.O 

IF < NTAU1 . EO • 1 >  GO  TO  85 

GO  TO  90 

DO  86  II  a  1 »NP 

DO  86  JJ  *  l »NWP 

PSTAUl < II fJJ>  a  o.O 

GO  TO  67 

DO  91  II  a  2rNTAUl 
12  a  NTAUi-Il+1 
DO  92  I  a  1 r NP 
DO  92  J  a  1 t NWP 

PSTAU1 <  I  r  J)  a  PSTAUl < I ? J )  +  PS(IfJ) 

DO  93  I  a  i,np 
DO  93  J  ~  1 f NWP 
PS(IfJ)  0.0 
DO  93  K  a  l f Np 

PS(I»J)  a  P8<I»J)  +  PHI T 1 < I r  K  r 1 2 ) *B1P ( K  f  J ) 
DO  94  I  =  1 f  NP 
DO  94  J  a  1 r NWP 
PS < I » J ) «DELHLF*PS < I f J ) 

DO  96  I  a  1 ,NP 

DO  96  J  If  NWP 

PSTAUl ( 1 ? J)  a  PS<IfJ)  +  PSTAUl < I f J) 

CONTINUE 

DO  68  I  a  ;i,np 

DO  68  J  a  l j NWP 

PS(IfJ)  a  o.O 

IF (NTAU.EO. 1 )  GO  TO  55 

GO  TO  69 

DO  58  II  a  i  ,NP 

DO  58  JJ  «  If NWP 

PSTAUIIIfJJ)  a  o.O 

GO  TO  77 
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69  DO  70  II  =  2i>NTAU 
12  a  NTAU--1 1  +  1 

DO  71  I  «  1 f NP 
DO  71  J  IfNUP 

71  PSTAU(IfJ)  =  PSTAU<I»J)  +  PS(ItJ) 

DO  72  I  =  1 f  NP 

DO  72  J  =  1»NWP 
PS(1»J)  =  0.0 
DO  72  K  a  1 r NP 

72  PS<I»J>  ~  PS(IrJ)  +  PHIT1 < I »K» I2)#B1P(K» J) 

DO  73  I  “  1 r NP 

DO  73  J  «  1 t NWP 

73  PS  <  I » J )  DELHLF#PS  <  I » J ) 

DO  76  I  =  1  *  NP 

DO  76  J  -  1 t NWP 

76  PSTAU<I»J)  «  PS(I»J)  +  PSTAU(IrJ) 

70  CONTINUE 

77  CONTINUE 

C 

C  ****************************************************************** 

C  START  TIME  LOOP 

C  TIME  »  0 

C  ***************************************************************** 

C 

NN1  a  1 

DO  412  I  a  If NOP 

412  UP(I)  =  W ( NN1 ) -YC1 < I ) 

DO  413  I  a  1 t NP 

Q  a  0.0 

DO  414  J  =  1 .NP 

414  Q  a  Q  +  PHTAUKI.  J)*XP1<  J) 

413  XWK.I)  =  Q 

DO  415  I  a  1 j  NP 

0  =  0.0 

DO  416  J  a  1 t NUP 

416  0  a  Q  +  PSTAUKI. J)*UP< J) 

415  XW2( I )  =  0 

DO  417  I  a  1 »NP 

417  XPTAU1  <  I )  a  XWKI)  +  XW2<  I ) 

DO  418  I  =  1 r NOP 

0  a  0.0 

DO  419  J  a  l,NP 

419  0  a  Q  +  CP(I» J)*XPTAU1< J) 

418  YPTAU1 < I )  a  Q 

DO  420  I  a  1 f NUP 

0  a  0.0 

DO  421  J  =  IfNC 

421  Q  =  Q  +  HC<I»J)*XC3<J> 

420  XWKI)  a  Q 

DO  422  I  a  1 » NUP 

Q  =  0.0 

DO  423  J  a  1 f NOP 
423  0  a  Q  +  EC< I » JXYPTAUl ( J) 

422  XW2<  Z )  =  Q 


1 

i 


424 


426 

425 


428 
427 

429 

4130 

C 

C 

C 


256 


257 


261 

260 

G 

t: 

c 

c 


301 

300 


303 

302 


304 


HO  424  I  =  1 f NUP 

YC3< I )  =  XWl(I)  +  XU2 ( 1 ) 

DO  425  I  =  1 f  NC 

0  *  0.0 

DO  426  J  »  1 r NC 
0=0+  FC(If J)#XC3<J) 

XWKI)  -  Q 
DO  427  I  ~  1 t NC 
Q  ~  0.0 

DO  428  J  ~  1 f NOP 
0=0+  GC<If J)*YPTAU1<J) 

XW2 ( 1 )  =  0 

DO  429  I  =  KNC 

XC3 ( I )  *  XWl(I)  +  XW2 ( I ) 

DO  430  I  =  1 » NUP 

E2 < I )  =  YC3( I )  ~  YC2 ( I ) 

DO  220  III  «  1» 300 

j*;******###******#************************#********************** 

CALCULATE  XP ( TIME+TAU ) f  AND  YP ( TIME+TAU ) 

tt*****************t******************************************** 

DO  253  I  =  1 t NP 

0  =  0.0 

DO  254  J  =  1 f  NP 
0  =  Q  +  PHTAUdf  J)*XP(  J) 

XWKI)  =  0 
DO  255  I  =  1 » NP 
0  -  0.0 

DO  256  J  =  1 j NUP 
0=0  +  PSTAUdf  J)#UP<J) 

XW2<  I )  ■-  0 

DO  257  I  »  1 f NP 

XPTAU(I)  =  XWld)  +  XW2 ( I ) 

DO  260  I  =  1» NOP 

0  =  0.0 

DO  261  J  =  1 f NP 
0=0+  CPdf  J)*XPTAU<  J) 

YPTAU(I)  =  0 

************************************.*****************■*********** 
TIME  =  TIME  +  T 

CALCULATE  XC2dIME+T+TAU)  f  AND  YC2 ( TIME+TAU) 

******JK***))C#****#!K*********#***)(C#******)(t*#*!(t*##*#**3((#*#**»*#**#* 

DO  300  I  =  1 f NUP 

0  =  0.0 

DO  301  J  ~  1»NC 
0=0+  HCdfJ)*XC2< J) 

XWKI)  =  0 
DO  302  I  =  KNUP 
0  ■  0.0 

DO  303  J  =  NOP 
0  =  Q  +  EC< I f J)#YPTAU< J) 

XW2( I )  =  Q 

DO  304  I  =  If NUP 

YC2< I )  =  XW1<I)  +  XW2( I ) 

DO  305  I  =  IfNC 
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( 

i 

1 

I 

«*> 

i 

V 

■j 


■j 

I 

I 

il 

1 

1 


a  =  o.o 

DO  306  J  =  IfNC 

306  Q  =  Q  +  FC<If J)*XC2< J) 

305  XWl(I)  =  Q 

DO  307  I  ■■=  1 1 NC 

a  ~  o.o 

DO  308  J  =  If  NOP 

308  a  «  a  +  GC(If J)#YPTAU< J) 

307  XW2 ( I )  =  0 

DO  309  I  =  IfNC 

309  XC2(I)  =  XUKI)  +  XU2< I ) 

C!*  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  j|^  ^  ^  j|^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  j|^  ^  ^  ^  ^  ^  ^  ^  ^  j|^  ^  ^  ^  .^.  ^  ^ 

C  CALCULATE  El 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  X  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  >L  ^  ^  ^  .  I .  ^  ^  ^  ^  ^  .K  ^  ^  J.  ^  J.  ^  ^  ^  ^  ^  ^  ■  i. 

^  ^  ^  ^  ^  ^  *  T*  T*  *  V  V  V  ^  *  T*  *  ^  ^  ^  ^  ^  <p  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

DO  290  I  »  1 t NUP 
290  El(I)  «  YC3  ( I )  -  YC2  ( I ) 

C  *************************************************************%*** 

C  ESTIMATE  XP*<TIME)f  AND  YP*(TIME> 

C,  *#*####*)«***###  **##*)f(****>|t###***)(!*)|t)(t#****#)|(*)tt)|C##*)K*))!)f!**)«  ###*##**# 

DO  432  I  =  If NUP 

432  YP2< I )  =  YP ( I )  -  YP1 < I ) 

DO  405  I  =  1 r  NP 

0  =  0.0 

DO  406  J  a  1 f NP 

406  Q  a  Q  +  PHITdf J)*XP1< J) 

405  XW1(I>  a  a 

DO  407  1  a  1 » NP 

a  =  o.o 

DO  408  J  a  l r NUP 

408  0  =  0  +  PSITdf  J)*UP<  J> 

407  XW2( I )  «  0 

DO  409  I  =  1 r NP 

0  =  0.0 

DO  410  J  =  1 » NUP 

410  0=0+  PSIT2(IfJ>*YP2< J) 

409  XW3( I )  a  0 

DO  411  I  a  i,NP 

411  XPKI)  ■  XUKI)  +  XW2<  I )  +  XW3  ( I  > 

DO  433  I  »  If NOP 

0  a  o.O 

DO  434  J  =  1  fNP 
434  0  =  0  +  CPdf  J)*XP1(J) 

433  YP1 < I )  =  0 

C  *#***#*#*#***#**#**####*###**#**###*##*#**#***##*#*#*****#*♦**## 

C  CALCULATE  XP(TIME)»  AND  YP(TIME) 

C  *#*****#**#*)(OK#*!«JK**###****#*#*#**##**###***#***)(t******J|t)([)|t##*)(t)(!# 

373  DO  500  I  a  1 » NP 

0  =  0.0 

DO  501  J  =  1  fNP 
501  0  «  Q  +  PHITdf J)*XP<J> 

500  XUKI)  =  0 

DO  502  I  -  1 f NP 

0  a  0.0 

DO  503  J  =  If NUP 
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503  Q  ~  Q  +  PSIT<If J)*UP< J) 

502  XW2<I)  *  Q 

DO  504  I  =  1 f  NP 

504  XP(I)  =  XWKI)  +  XU2( I ) 

DO  507  I  =  1 r NOP 

Q  *  0.0 

DO  508  J  =  1 f  NP 
508  Q  “  Q  +  CP(IrJ)*XP(J) 

507  YP  < I )  =  Q 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  X  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  J.  ^  ^  ,  4«-  ^  ^  ^  .1.  ^  ^  ^  4^  ^  ^  ^  4^  4k  ^  ^  4^  4^  ^  ^  ^  ,1. 

^  ^  ^  4^  4^  ^  4^ 

C  CALCULATE  XCl<2*TIME)f  AND  YCl(TIME) 

C  ####*#!)(**  Jtt####*#*###**####***##*#####**#*###*####**##*##*##**#** 

DO  700  I  *  IfNUP 

Q  =  0.0 

DO  701  J  *  1 f NC 

701  Q  “  Q  +  HC ( I r  J ) #XC1 <  J  > 

700  XWKI)  =  Q 

DO  702  I  *  1 f  NUP 

a  -  o.o 

DO  703  J  =  If NOP 

703  Q  *  0  +  EC< I f J)#YP( J) 

702  XW2 ( I )  ~  Q 

DO  704  I  =  IfNUP 

704  YCl(I)  "  XWKI)  +  XW2  ( I ) 

DO  705  I  -•  1  f  NC 

Q  =  0,0 

DO  706  J  ~  IfNC 

706  Q  =  Q  +  FC(If J)*XC1(J) 

705  XWKI)  ~  Q 

DO  707  I  •  KNC 

Q  »  0.0 

DO  708  J  -  If NOP 

708  Q  =  Q  +  GC( I f J)#YP< J) 

707  XW2 ( I )  *  0 

DO  709  I  ~  IfNC 

709  XCl(I)  ~  XWKI)  +  XW2  ( I ) 

( j  )|C  JjC  JjC  *  3^(  ]|[  JfC  jjC  )(t  j|(  ]|c  jf[  )jc  3jt  ]|c  ]|(  ])( $  j|c  ]|C  3J(  )|(  ]|(  j|C  jf(  ])( 5jC  ]|(  )|E  )|(  5)£  ]|{  )|(  J(C  3fC  )|E  Jj( )(( )(!  3jC  )|E  ]|C  JjC  IfE  )(C  I|t  Jjt  JjC  ]|C  Jjt  E|[  <(t  E(E  ))C  3j(  J|E  EjE  3)t  I^E  ]|C  Ejt 

C  ESTIMATE  XP* ( TIME+TAU ) f  AND  YP*(TIME+TAU) 

NN1  =  NN1  +  1 
DO  221  I  »  If  NOP 
221  UP(I)  =  W(NN1 )-YCl < I ) 

DO  520  I  ■  IfNP 

Q  =  0.0 

DO  521  J  =  IfNP 

521  Q  “  Q  +  PHTAUKIf J)*XP1<J) 

520  XWKI)  =  Q 

DO  522  I  =  1»NP 

Q  -  0,0 

DO  523  J  •---  IfNUP 

523  0  -  Q  +  PSTAU1 < I r J) #UP ( J) 

522  XW2< I )  =  Q 

DO  524  I  .»  IfNP 

524  XPTAUKI)  -  XWKI)  +  XW2<I> 
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1  .  NOP 


DO  S25  I  = 

Q  =  0.0 

DO  526  J  =  1  .  NF' 

526  Q  ~  Q  +  CP<I.J)*XPTAU1< J) 

525  YPTAU1 ( I )  “  Q 

A  A  ^  ^  U/  ^  U/  X  ^  ^  ^  dr  dr  d  u'  X  ^  |L  «d  U/  ^  d-  dr  d/  d#  df  d>  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d 

^  ^  ®  ^  ^  ®  ^  ^  ^  ^  ®  ®  ^  ®  ^  ®  ^  ^  t*  v  T  ^  ®  ^  ^  ^  ®  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^p  ^  ^ 

C  CALCULATE.  XC*(2*TIHE+TAU*) .  AND  YOMTIME+TAU*) 

A  A  A  4[  ^  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  ^  ^  d  4  d  d  d  d  d  d 

*  ^  ^  ^  V  ^  ^  V  ^  ^  ^  ^  ^  ^  *  *  *  V  T*  *  V  v  O’  fl*  *v  *  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^p  ^  ^  ^  ^p  ^  ^  ^ 

DO  529  I  =  1 . NUP 

Cl  =  0.0 

DO  530  J  -=••  J.  »  NC 

530  Q  ~  Q  +  HC  < I . J ) *XC3 ( J ) 

529  XW1 < I )  =  Q 

DO  531  I  «  1 » NUP 
Q  a  0.0 

DO  532  J  “=  1  .  NOP 

532  Q  =  Q  +  EC  < 1 . J ) #YPTAU1 ( J ) 

531  XW2 ( I )  =  Q 

DO  533  I  a  1 . NUP 

533  YC3 ( I )  a  XWI(I)  +  XU2 < I ) 

DO  535  I  =  1 » NC 

Q  =  0.0 

DO  536  J  »  l.NC 
536  Q  =  Q  +  FC<If J>*XC3< J) 

535  XWKI)  =  Q 

DO  537  I  a  1 » NC 
Q  =  0.0 

DO  538  J  a  1 . NOP 

538  Q  =  Q  +  GC<I.J)*YPTAU1<J) 

DO  539  I  a  l.NC 

539  XC3  ( I )  =  XWKI)  +  XW2(I) 

C  CALCULATE  E2 

Q  ^  ^  ^  ^  ^  ^  ^  y|^  ^  ^  ^  ^  y^p  ^  y^  ^  yj^  y^  ^  ^  ^  y^p  ^  y^p  y^p  y|^  ^  y^  ^  ^  ^||  y^  y^  y^  y^p  ^  yj^  y^  y^p  y^  y^r  y^r  y^p  y|^  y|^  ^  j^p  y^p  j|^  j|p  y^  y^  p^p 

DO  540  I  a  1 .  NiJP 

540  E2 ( I )  a  YC3 ( I )  -  YC2  < I ) 

DO  560  I  =  1  .  NUP 
EA(IIl)  =  El(I) 

560  EB(IIl)  a  E2( I > 

220  CONTINUE 

SMEANA  a  o.O 

SMEANB  a  o.O 

DO  561  I  =  201.300 

SMEANA  a  SMEANA  +  E.A<I) 

561  SMEANB  =  SMEANB  »  EB(I) 

SMEANA  =  SMEANA/100 

SMEANB  a  SMEANB/100 

VEA  =  o.O 

VEB  =  o.O 

DO  563  I  a  201.300 
VEA  a  VEA  +  <EA< I) -SMEANA) **2 
563  VEB  *  veB  +  < EB< I ) -SMEANB > **2 
PEASS<KK2)  a  VEA/100 
PEBSS(KK2)  a  VEB/ 100 
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0--  <J1 


1800  CONTINUE 

570  FORMAT ( 5X » 'PEASS<  TAU*  =  ',F12.8 
WRITE<6  *  571 ) < PEASS<  I )  r 1-1 » 6 ) 

571  FORMAT (5X»F18.10) 

WRITE <6, 572)  TAU1 

572  FORMAT <5Xr'PEBSS<  TAU*  =  ',F12.8 
WRITE  <  6  f 571 ) (PEBSS ( I > » I~1 » 6 ) 

1801  CONTINUE 
STOP 
END 

SUBROUTINE  RANDU ( IX > I Y , YFL ) 

IY  =  IX*6553? 

IF<IY>5»6,6 

IY  ~  IY  +  2147483647+1 
YFL  *  IY 

YFL  =  YFL#0«  46566 13E-9 
RETURN 
END 
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EXAMPLE - 2TH  ORDER  PLANT.  1ST  ORDER  CONTROLLER 

NO.  OF  PLANT  STATES  =  2 
NO.  OF  PLANT  INPUTS  =  1 
NO.  OF  EXTERNAL  INPUT  =  1 
NO.  OF  PLANT  OUTPUTS  =  1 

NO.  OF  CONTROLLER  STATES  (  EACH  CONTROLLER)  =  1 

PLANT  STATE  MATRIX  —  AP 

.0  1.0 

.0  -10.0 

PLANT  CONTROL  INPUT  MATRIX  —  BP 

.0 

200.0 


PLANT  OUTPUT  MATRIX  —  CP 

1.0  .0 


CONTROLLER  STATE  MATRIX  —  FC 
.523810 


OBSERVER  MATRIX  —  GE 
20 . 0 
25.0 

CONTROLLER  CONTROL  INPUT  MATRIX  —  GC 
-.18162 

CONTROLLER  OUTPUT  MATRIX  (STATES)  —  HC 

1.0 


CONTROLLER  OUTPUT  MATRIX  (INPUTS)  —  EC 
1.381 
NT  *  51 

T  -  SAMPLE  PERIOD  *  0.0125  SEC 

DELTA  *  T/(NT~1 )  -  INCREMENT  USED  IN  THE  NUMERICAL 

INTEGRATIONS  TO  COMPUTE  PSITAU.PSI T .PSIT 
PSITAU1  USING  TRAPEZOIDAL  RULE. 

U  IS  THE  EXTERNAL  INPUT  (WHITE  GAUSSIAN  NOISE  WITH 
MEAN  =0.0.  AND  VARIANCE  «  1.0 
THE  STEADY  STATE  SAMPLE  VARIANCE  OF  ERRORS 
PEASS(  TAU*  «  0.0  ) 

0.0 

0.0002195683 
0.0008547062 
0.0018880414 
0.0033245913 
0.0051899776 
PEBSS (  TAU*  =  0.0  ) 

0.0051899 776 
0.0033477221 
0.0019137899 
0.0008709673 
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0.0002244494 

0.0 


PEASS (  TAU*  =  0.0025  ) 
0.0002195683 
0.0 

0.0002097336 
0.0008310268 
0.0018685847 
0.0033477221 
PEBSS (  TAU*  =  0.0025  ) 
0.0074686036 
0.0052524656 
0.0034256792 
0.0019714807 
0.0008955190 
0.0002240368 

PEASS <  TAU*  »  0.005  > 
0.0008547062 
0.0002097336 
0.0 

0.0002076342 
0.0008371675 
0.0019137899 
PEBSS <  TAU*  *  0.005  > 
0.0100931898 
0.0075152330 
0.0053092465 
0.0034589050 
0.0019702637 
0.0008699684 

PEASS (  TAU*  -  0.0075  ) 
0.0018880414 
0.0008310268 
0.0002076342 
0.0 

0.0002127137 
0,0008709673 
PEBSS (  TAU*  -  0.0075  > 
0.0130441934 
0.0101163760 
0.0075445175 
0.0053128712 
0.0034278994 
0.0019165913 

PEASS (  TAU*  =0,01  ) 
0.0033245913 
0.0018685847 
0.0008371675 
0.0002127137 
0.0 
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0.0002244494 

PEBSS<  TAU*  =  0.01  ) 
0.0163251571 
0.0130585469 
0.0101340003 
0.0075355470 
0.0052701645 
0.0033652126 

PEASS<  TAU*  *  0.0125  ) 
0.0051899776 
0.0033477221 
0.0019137899 
0.0008709673 
0.0002244494 
0.0 

PEBSS<  TAU*  ~  0.0125  ) 
0.0199597441 
0.0163656212 
0.0131007358 
0.0101497732 
0.0075195357 
0.0052378476 


APPENDIX  F 


COMPUTER  PROGRAM  LISTING 
FOR 

ALGORITHM  FOR  ESTIMATING  t 
AND  EXAMPLE  OF  OUTPUT  WRITTEN 
IN  FORTRAN 
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DIMENSION  AP(2*2)»B1P(2>2)  ,GE( 2»2)  >  CP ( 1  > 2  >  ,  FC  (2»  2 )  >GC <  2 » 1 )  r 

1  HC(2*2) ,EC(2rl> »PHITl(4f4x 101) >PSIT1(4, 4) rPHTAU(4,4) ,YP2(2> r 

2  PHIT ( 4* 4 )  > PSTAU( 4»4) »PSIT(4>4) r 

4  INDEX (4) » W ( 4000) ,PS ( 4 , 4 ) , PHTAU1 ( 4 , 4 ) » PSTAU1 ( 4 , 4 > , YC1 ( 2 ) » 

5  YC2<2>  rEl(2)>E2(2) ,XP(2) rXFTAU(2> »XP1(2> r YP < 2 > » YPTAU ( 2 >  * 

6  YP1 ( 2 ) » XC1 ( 2 ) »  XC2<  2)  »AM(4»4)*F’T(4»4)»Pl(4r4)rDl(4»4)  >■ 

7  D2( 4 » 4 ) » D3( 4 )  »XW1  (2)  »XW2(2>  »ECCP(2»  2)  »D(4f4)»  YPI AIJ1  (2 )  r 
S  PEASS(SO) fPEBSS(SO) »EA(300) »EB(300) > YC3(2> »XC3(2) »UP(2) v 
9  XW3<2) rXPTAUl (2) 

CCCCC  PROVIDE  MAXIMA  FOR  CALLED  ARRAYS 
NPM  *  2 
NUPM  ~  2 
NWPM  *  2 
NOPM  =  1 
NCM  »  2 

NHM  ■  NPM  +  NUPM 

NFM  *  NPM  +  2*NCM 

NRRM  »  2*NPM  +4 

C 

^  ^  ^  ^  X  ^  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  L  X  x  X  X  X  X  X  X  J  'I  ■  I 

^  ^  ^  ^  ®  ^  ^  ^  ^  ^  ®  ^  ^  ^  ^  ^  ^  ^  ^  ^  rp  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  /|\  ^ 

C  READ  INPUT  DATA 

C 

WRITE (6, 899) 

899  FORMAT ('1') 

100  REAti (  5 » 900 )  ID 
900  FORMAT (20A 4) 

WRITE (6» 902)  ID 
902  FORMAT  (  '  1 '  >  20A4 ) 

READ ( 5  r 906  >  NP » NUP  » NWP . NOP  r  NC 
906  FORMAT (513) 

WRITE ( 6 , 908 )  NP - NUP  r NWP  *  NOP » NC 
908  FORMAT  ( 'ONO.  OF  PLANT  STATES  *  'j-13/ 

1  '  NO.  OF  PLANT  INPUTS  *  '  *13/ 

2  '  NO.  OF  DISTURBANCE  INPUTS  «  ' »  13/ 

4  '  NO,  OF  PLANT  OUTPUTS  «  ' ,  13/ 

5  '  NO.  OF  CONTROLLER  STATES  (EACH  CONTROLLER)  «  '*13) 

WRITE( A* 910) 

910  FORMAT ( 'OPLANT  STATE  MATRIX  —  AP'> 

110  DO  112  I  =  1 »NP 

READ( 5»914  >  (AP(I»J)»J*1*NP) 

112  WRITE(A»913>  (AP(I»J)*J=1»NP> 

913  FORMAT ( '  '*8G13.A) 

914  FORMAT ( AFl 2 . 7 ) 

915  FORMAT (861 3. A) 

WRITE( A»91A> 

916  FORMAT ( 'OPLANT  CONTROL  INPUT  MATRIX  —  B1P') 

120  DO  122  I  *  1 »NP 

READ(5»914)(B1P(I»J)*J=1* NUP ) 

122  WRITE(A*913XBiP(I,J>,J*i,NUP) 

WRITE(A*918) 

918  FORMAT <  'OOBSERVER  MATRIX—  6E') 

130  DO  132  1*1 »NP 

READ(5»914) (0E( I » J) » J*1 »NWP) 
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132  WRITE(6f913)(GE(IfJ)fJ=1 fNWF) 

WRITE(6f920> 

920  FORMAT ( 70PLANT  OUTPUT  MATRIX  —  CP') 

140  DO  142  1=1 fNOP 

READ(5f914) (CP( I f J) f  J=1 fNP ) 

142  WRITE(6f913) (CP( I f J) f J=1 fNP) 

WRITE(6f922) 

922  FORMAT (  "OCONTROLLER  STATE  MATRIX  ~  FC7> 

150  DO  152  I  =1fNC 

READ(5f914XFC(IfJ)  f  J=1fNC> 

152  WRITE(6f913XFC(IfJ)fJ=1fNC) 

WRITE(6f924) 

924  FORMAT < 'OCONTROLLER  CONTROL  INPUT  MATRIX  ~  GC  ') 

160  DO  162  1=1 fNC 

READ (5f914) (GC(IfJ>fJ=1 fNOP) 

162  WRITE  (6  f  913)  (GC  <IfJ)fJ=!1  fNOP  ) 

WRITE(6f925) 

925  FORMAT < 7 ©CONTROLLER  OUTPUT  MATRIX  (STATES)  ~  HC7) 

170  DO  172  1=1fNUP 

READ (5f 914) <HC(IfJ)fJ=1fNC) 

172  WRITE(6f913XHC(IfJ)fJ=1fNC> 

WRITE(6f92 6) 

926  FORMAT < 7 OCONTROLLER  OUTPUT  MATRIX  (INPUTS)  —  EC7) 

180  DO  182  I=1fNUP 

READ(5f914)(EC(IfJ)fJ=1fN0P) 

182  WRITE(6»913XEC(IfJ).J=1»N0P> 

READ(5f928)  TfNT 
928  FORMAT (F10.4fI5) 

XNT  =  NT 

DELTA  =  T/(XNT-1) 

WRITE(6f930)  TfNT 

930  FORMAT ( 7  IT  ■  7fF10.4/ 

1  7  NT  ■  7  f 15/ 

3  7  T  ■  SAMPLE  RATE,7/ 

4  7  NT-1  ■  NB.  OF  EVENLY-SPACED  SUBINTERVALS  INTO  WHICH  T  IS'/ 

5  7  DIVIDED.7/ 

6  7  DELTA  =  T/<NT-1 )  «  INCREMENT  USED  IN  THE  NUMERICAL'/ 

7  7  INTEGRATIONS  TO  COMPUTE  VZTr  VZTAUf  AND7/ 

9  7  VZTAU1  USING  THE  TRAPEZOIDAL  RULE.7//) 

WRITE(6f931 ) 

931  FORMAT <2Xf  7  W  IS  THE  DISTURBANCE  VECTOR  (WHITE  GAUSSIAN  NOISE7/ 

1  7  WITH  MEAN  ■  0  AND  VARIANCE  *  l)7) 

C 

C  ***************************************************************** 

C  GENERATE  WHITE  GAUSSIAN  NOISE  WITH  MEAN  *  0  AND  VARIANCE  *  1 
C  ***************************************************************** 

c 

ix  *  mil 

DO  192  I  *  If  1200 
A  »  0*0 

DO  193  J  «  1 f 12 
CALL  RANDU(IXfIYfY) 

IX  »  IT 
193  A  -  A4Y 
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192  W(I)  -  A-6 

DO  1199  I  =  DNUP 
DO  1199  J  *  1,NP 
ECCP(IrJ)  «  0.0 
DO  1199  K  =  1 t NOP 

1199  ECCP(IfJ)  *  ECCP < I f J )  +  EC(IrK)»CP<K'J> 

802  FORMAT (5X» 50 13. 6) 

835  FORMAT <5X> 5013. 6) 

C 

C  **#*####****##*#####*##***#***#*****#***###****♦#*#*****:*******.** 

C  CALCULATE  PHIT(T)»  PSIT(TAU>» 

C  PH  IT  (0)  r  PHIT  (DELTA)  »  PHITt  2DELTA)  , .  ♦  .  »PHIT'T  DEL  I A )  » 

C  PS IT < 0 > »  PS IT (DELTA) »  PS IT < 2DELTA )»♦..» PSIT ( T-  DELTA > 

C  ****************$******$$#$*****$**Ht*******>t$*#**t*#*t*t*  #***){*  .*.* 

C 

TI  a  0.0 

DELHLF  a  DELTA/2.0 
DO  402  II  *  1>NP 
DO  402  JJ  a  1 r NP 
IF(II.EQ.JJ)  GO  TO  403 
PHIT1 ( 1 1 f  JJ » 1 )  a  0.0 
GO  TO  402 

403  PHITKII,  JJrl)  =  1.0 

402  CONTINUE 

DO  4  II  a  2  r  NT 
TI  a  TH DELTA 
PHITKblrlt)  a  i. 

PH I T 1 < 1 » 2  > 1 1 >  =  (1./10. )*(1.-EXP(-10.*TI)> 

PHIT1 <  2  >  1 » 1 1  >  »  0.0 
4  PHITl(2f2f  ID  =  EXP(-10.*TI) 

DO  400  II  =  1 r NP 
DO  400  JJ  a  i„NP 

400  PH I T < 1 1 » J J >  a  PHIT1 ( 1 1 » JJ > NT ) 

NTAU2  a  o 
N7  a  2000 
N8  =  300 
N9  a  50 
I COUNT  a  o 
I CUNT 1  a  1 
N4  a  0 
N5  a  51 
NTAU  «  4 
N3  a  <N4+N5)/2 
600  NTAIJ1  a  N3 

I COUNT  a  I COUNT  +  1 

CALL  COVAR<AP»BlP»CP>FC»GC»HCrECrGE»ECCP» WrPHIT »PHIT1 f NTAU r NTAIJ 1 
1  EA » PEASS  f NPM » NUPM  t NOPM  t NCM » NHM » NWPM » N7  r  N8 » I COUNT » N9 » DELHLF ) 

WRITE <6» 655)  NTAUDNTAU2 
655  F0RMAT<5XfI5>5X>I5) 

NERROR  a  IABS ( NTAU1-NTAU2 ) 

IF(NERROR«EQ. 1 )  GO  TO  605 
IF < ICUNT1 .NE . 1 )  GO  TO  601 
N3  =  N3  +  1 
NTAU3  a  N3 
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NTAU1  =  NTAU3 
ICUNT2  =  ICOUNT 
ICUNT 1  =  0 
GO  TO  600 

601  IF(PEASS( ICUNT2) *LT *PEASS( ICOUNT) >  GO  TO  602 
N4  =  N3  -  1 

N3  ««  <  N4  +  N5  )/2 
NTAU2  =  NTAUi 
ICUNT1  =  1 
GO  TO  600 

602  N5  -  N3  -  1 

N3  -  (N4+N5 )/2 
NTAU2  =  NTAUI 
I CUNT 1  =  1 
GO  TO  600 

605  IF (NTAUI .  EQ. 1 >  GO  TO  606 

IF (NTAUI .  LT . NTAU2)  GO  TO  608 
ICUNT2  =  ICOUNT  -  1 

IF(PEASS( ICUNT2) .LT .PEASS ( ICOUNT) )  GO  TO  613 
616  NTAUI  ■  NTAUI  +  1 

ICOUNT  »  ICOUNT  +  1 

CALL  COVAR( AP » B1P*CP»FC»  GC»HC  *EC»GE»ECCP» W*PHIT  *PHIT1 r NT AUr NTAUI 
1  EA  * PEASS » NPM * NUPM , NOPM * NCM *  NHh » NWPM  *  N7  *  N8 * I COUNT * N9  r  DELHLF > 
ICUNT2  =  ICOUNT  -  1 

IF ( PEASS (ICUNT2>*LT. PEASS (ICOUNT))  GO  TO  615 
TAU1  =  (NTAU 1-2) 40 *0125/50 
GO  TO  612 

615  IF(NTAU1.LT.N5>  GO  TO  616 
TAU1  *  (NTAU1-1)*0. 0125/50 
GO  TO  612 

613  TAU1  =  (NTAUl-2>*0. 0125/50 
GO  TO  612 

606  ICUNT2  =  ICOUNT  -  1 

IF ( PEASS ( ICUNT2) *LT • PEASS ( ICOUNT ) )  GO  TO  607 
TAU1  =  ( NTAUI -1 ) #0*0125/50 
GO  TO  612 

607  TAU1  =  NTAU140. 0125/50 
GO  TO  612 

608  ICUNT2  =  ICOUNT  -  1 

IF ( PEASS ( ICUNT2) *LT • PEASS ( ICOUNT ) )  GO  TO  609 
NTAUI  «  NTAUI  -  1 
ICOUNT  «  ICOUNT  +  1 

CALL  COVAR( AP*B1P»CP*FC  *GC*HC*EC*GE*ECCP» W*PHIT rPHITl *NTAU* NTAUI 
1  EA* PEASS  *  NPM  »  NUPM  *  NOPM  *  NCM  *  NHM »  NWPM  *  N7  *  N8  r ICOUNT » N? *  DELHLF ) 
ICUNT2  »  ICOUNT  -  1 

IF ( PEASS <ICUNT2)*LT. PEASS (ICOUNT))  GO  TO  610 
TAU1  -  NTAU1B0. 0125/SO 
GO  TO  612 

610  TAU1  *  (NTAU1-1 >*0.0125/50 
00  TO  612 

609  TAU1  -  ( NTAUI -1 >40.0125/50 

612  WRITE(6*650>  TAU1 

650  FORMAT (5Xf '  THE  SKEW  OF  THE  SECOND  CONTROLLER  IS  ■  '»F13.8> 

657  FORMAT <  5X  » ' THE  NUMBER  OF  ITERATIONS  18  '*I5> 
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WRITE<6»657>  ICOUNT 

STOP 

END 

SUBROUTINE  COVAR <  AP  f  B1P  »  CP » PC » GC » HC » EC • GE r ECCP »  U » PHI  I .PHIT1. 

1  NTAUn NTAU1 » EA » PEASS , NPM » NUPM * NOPH » NCM » NHH » NWPM  » N7  >  N8  r ICOUNT >  N9 

2  DELHLF ) 

DIMENSION  ECCP < NPM » NUPM) 

DIMENSION  AP(NPM»NPM> »B1P<NPM»NUPM> *CP(NOPM»NPM) »FC(NCM»NCM> > 

1  GC ( NCM » NORM )  » HC  < NUPM » NCM ) » EC <  NUPM » NOPM )  »GE< NPM  y  NWPM ) , W  <  N7 ) • 

2  PHIT <NHM»NHM )  »PHIT1  ( NHM»NHM » 101  >  »EA(NS>  >F‘EASS < N9 > » 

3  PHTAU( 4  » 4  >  »PHTAU1 ( 4 » 4) xPSITl <  484  >  >  PSTAU (4*4) fYCl<2) *  YP2  <  2 ) t 

5  YC2<2>  *E1 (2) »E2(2>  »XP(2) »XPTAU(2) *XP1<2> »YP<2> »YPTAU(2) * 

6  YP1(2)*XC1<2) *XC2 (2) * AM( 4 * 4 ) n PS < 4 * 4 ) r PI ( 4 * 4 ) »D1 ( 4 *4 ) , 

7  D2 ( 4  *  4 )  *  D3 ( 4 )  *XW1  (2)  *XW2<  2 )  *  D<  4  *  4 )  r  YF'TAUl  ( 2) » YC3  ( 2) r  XC3 ( 2 >  * 

9  XW3(2> > XPTAU1 (2) »PSTAU1 ( 4 * 4 ) ,PSIT2 ( 4 , 4 ) *PSIT<4* 4 ) » UP ( 2) 

NF  «  2 
NOP  a  1 
NC  «  1 
NUP  =  1 
NUP  a  1 
NT  a  50 

DO  31  I  a  1 rNUP 
YC2 ( I )  a  o.o 
YC3 ( I )  a  O.o 

31  YC1 < I )  =  0.0 

DO  32  I  a  l,NUP 

32  El (I)  a  YC3< I )  -  YG2< I ) 

C  FROM  INITIAL  VALUE  XP(TIME)r  AND  YP ( TIME )  ARE  EQUAL  TO  ZERO 

DO  35  I  a  1 , NP 
XPl(I)  a  0,5 

35  XP(I)  a  0.0 

DO  435  I  a  1 f NOP 
«  a  0.0 

DO  436  J  a  1 f NP 
436  0  =  Q  +  CP<I?J)*XP1< J) 

435  YPl(I)  a  Q 

DO  36  I  =•  1»  NOP 

36  YP ( I )  ~  0.0 

C  FROM  INITIAL  VALUE  XCl < TIME+T ) *  YC1 ( TIME )  ARE  EQUAL  TO  ZERO 
DO  37  I  a  l f NC 
XCl < I )  a  0,0 
XC3( I )  a  0.0 

37  XC2( I )  a  o.O 

DO  52  I  =  l,NUP 
52  YCl(I)  =  0.0 

IDEL  =  0 

DO  5000  I  =  i,nt 
IDEL  a  i DEL  +  1 

IF(NTAU.EQ.NTAU1 . AND. IDEL.EQ.NTAU)  00  TO  7 
IF( IDEL.EQ.NTAU)  GO  TO  16 
IF<ID€L.E0.NTAU1>  GO  TO  22 
GO  TO  5000 
7  DO  8  II  =  l,NP 

DO  8  JJ  «  i,NP 
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PHTAU<  II f JJ)  =  PHITKIIfJJfIDEL) 

8  PHTAUKIIf  JJ)  »  PHITK  IIf  JJf  I  DEL) 

GO  TO  5000 

16  DO  17  II  =  IfNP 
DO  17  JJ  =  IfNP 

17  PHTAU(IIfJJ)  =  PHITKIIfJJfIDEL) 

GO  TO  5000 

22  DO  23  II  =  IfNP 
DO  23  JJ  ~  IfNP 

23  PHTAUKIIf  JJ)  »  PHI  T1  <  1 1  f  J  J  f  I  DEL  > 

5000  CONTINUE 

DO  550  I  *  IfNP 
DO  550  J  =  IfNWP 
PKIfJ)  =  0.0 

550  PS(IfJ)  =  0.0 
DO  551  I  »  IfNP 
DO  551  J  =  1 r NWP 
PSIT2( I f  J)  a  0.0 
PSTAUCIfJ)  a  0.0 
PSTAU1 <  I  f  J)  =  0.0 
PSIT(IfJ)  a  0.0 
DO  551  K  =  IfNP 

PSIT2< I » J)  ■  PSIT2< I » J)  +  PHIT1<IfKfNT)*GE<NfJ) 
PSTAU(IfJ)  =  PSTAU(IfJ)  +  PHIT1<IfKfNTAU)#B1P<KfJ) 
PSTAU1  <  I  f  J)  =  PSTAUKIfJ)  +  PHIT1  ( I  fKfNTAUI  )*B1P(Kf  J) 

551  PSIT(IfJ)  =  PSIT(IfJ)  +  PHIT1 < I ,KfNT ) *B1P<Kf J) 

DO  552  I  »  1 fNP 

DO  552  J  =  IfNWP 

PSIT2< I f J)  a  DELHLF*PSIT2(If J) 

PSTAU(IfJ)  =  DELHLF*PST  AU  < I » J ) 

PSTAUKIfJ)  =  DELHLF*PSTAUK  I  f  J) 

552  PSIT(IfJ)  =  DELHLF*PSIT<If J) 

60  DO  61  II  =  2 f NT 
12  =  NT-I1+1 

DO  62  I  =  IfNP 
DO  62  J  =  IfNWP 

PSIT2<  I  f  J)  a  PSIT2<  I  f  J)  +  PKIfJ) 

62  PSIT(IfJ)  a  PSIT(IfJ)  +  PS<IfJ) 

DO  63  I  a  1,NP 

DO  63  J  =  IfNWP 
PKIfJ)  =  0.0 
PS(IfJ)  a  o.O 
DO  63  K  •  IfNP 

PKIfJ)  «  Pl<IfJ)  +  PHIT1  <  I  fKf  12 )  *GE  <Kf  J) 

63  PS(IfJ)  =  PS(IfJ)  +  PHITK I fKf I2)*B1P<K f J) 

DO  64  I  a  l,NP 
DO  64  J  a  IfNWP 

PKIfJ)  *  DELHLF*P1(IfJ) 

64  PS(IfJ)  *  DELHLF*PS<IfJ) 

DO  66  I  »  IfNP 
DO  66  J  a  IfNWP 

PSIT2(IfJ)  »  PKI,J)  +  PSIT2<IfJ) 

66  P8IT(IfJ)  *  PS ( I f J )  +  PSIT(IfJ) 

61  CONTINUE 
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DO  80  I  =  IfNP 
DO  80  J  =  IfNWP 
PS(IfJ)  -  0,0 
IFINTAUi.EG.l)  GO  TO  85 
GO  TO  90 
DO  86  II  a  1 r NP 
DO  86  JJ  ~  1 f NWP 
PSTAU1 (II » JJ)  ~  0.0 
GO  TO  67 

DO  91  II  a  2fNTAU1 
12  -  NTAU1-I1+1 
DO  92  I  =  1 » NP 
DO  92  J  ~  IfNWP 

PSTAUKIfJ)  =  PSTAUKIfJ)  +  PS(IfJ) 

DO  93  I  =  1 , NP 
DO  93  J  =  IfNWP 
PS<IfJ)  =  0.0 
DO  93  K  »  1 f  NP 

PS(IfJ)  -  PS(IfJ)  +  PHIT1<IfKfI2)*B1P(Kf J) 
DO  94  I  a  1„NP 
DO  94  J  ~  IfNWP 
PS  < I f  J ) aDELHLF*PS ( I f  J ) 

DO  96  I  a  l t NP 
DO  96  J  a  1,NWP 

PSTAUKIfJ)  a  PS(IfJ)  +  PSTAUKIfJ) 

CONTINUE 

DO  68  I  a  l fNP 

DO  68  J  a  l,NWP 

PS  <  I  F  J )  a  0.0 

IF (NTAU.EQ. 1 )  GO  TO  55 

GO  TO  69 

DO  58  II  a  IfNP 

DO  58  JJ  a  l,NWP 

PSTAU(IIfJJ)  a  0.0 

GO  TO  77 

DO  70  II  a  2  f  NT  All 
12  a  NTAU-Il+1 
DO  71  I  a  IfNP 
DO  71  J  a  i,NWP 

PSTAU(IfJ)  a  PSTAU(IfJ)  +  PS(IfJ) 

DO  72  I  -  IfNP 
DO  72  J  a  IfNWP 
PS(IfJ)  a  o.O 
DO  72  K  a  IfNP 

PS(IfJ)  =  PS(IfJ)  +  PHIT1(IfKfI2)*B1P(Kf J) 

DO  73  I  a  i,NP 
DO  73  J  =  IfNWP 

PS(IfJ)  =  DELHLF*PS<IfJ> 

DO  76  I  a  l,Np 

DO  76  J  a  l,NWP 

P5TAU< I f J)  a  PS(I,J)  +  PSTAU(IfJ) 

CONTINUE 
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C  START  TIHE  LOOP 

TIME  =  0 

C 

77  NN1  =  1 

DO  412  I  =  1 t NOP 

412  UP(I)  =  W(NN1 )-YCl ( I ) 

DO  413  I  »  1 f  NP 

Q  =  0.0 

DO  414  J  =  1 f  NP 

414  0  =  0  +  PHTAU1 < I » J ) #XP1 <  J ) 

413  XWl(I)  =  Q 

DO  415  I  =  1 f  NP 

Q  ■  0.0 

DO  416  J  »  IfNUP 

416  0  =  0  +  PSTAU1 ( I » J)#UP( J> 

415  XU2  < I )  =  0 

DO  417  I  =  1 t NP 

417  XPTAUl(I)  =  XU1<I)  +  XW2(I) 
DO  418  I  =  1 r  NOP 

0  ■  0.0 

DO  419  J  =  1 f  NP 

419  0=0+  CP< I f J)*XPTAU1 < J) 

418  YPTAU1 < I )  «  0 

DO  420  I  =  IfNUP 
0  =  0.0 

DO  421  J  ~  IfNC 

421  0=0+  HC<I» J)*XC3< J) 

420  XWl(I)  =  0 

DO  422  I  =  1 t NUP 

0  =  0.0 

DO  423  J  =  1 »NOP 

423  0=0+  EC(Ir J)*YPTAU1< J) 

422  XW2< I )  =  0 

DO  424  I  =  IfNUP 

424  YC3< I )  =  XW1 < I )  +  XW2 ( I > 

DO  425  I  =  IfNC 

Q  =  0.0 

DO  426  J  =  IfNC 

426  0  ■  Q  +  FC<If J)*XC3U> 

425  XW1 < 1 )  =  0 

DO  427  I  =  IfNC 

0  -  0.0 

DO  428  J  =  If NOP 

428  Q  ■  Q  +  GC ( I f  J ) *YPT AU 1 ( J ) 

427  XW2( I )  =  0 

DO  429  I  =  1»NC 

429  XC3<  I )  -  XWKI)  +  XU2  '  I ) 

DO  430  I  -  IfNUP 

430  E2< I >  *  YC3 ( I >  -  YC2<I) 

DO  220  III  =  1 f 60 

DO  253  I  =  1 f NP 

Q  *  0.0 

DO  254  J  *  1  f  NP 
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254  Q  =  Q  +  PHTAU<If J)#XP( J) 

253  XWKI)  =  G 

DO  255  I  ~  l.NP 

Q  ■  0.0 

DO  256  J  ~  1 t NUP 

256  G  *  Q  +  PSTALKIt J)*UP<J) 

255  XW2< I )  a  Q 

DO  257  I  -  l.NP 

257  XPTAU(I)  •»  XWl(I)  +  XW2 ( I ) 
DO  260  I  »  UNOP 

G  -  0.0 

DO  261  J  ~  1 »NP 
261  G  ~  G  +  CP<  I  f J)*XPTAU< J) 
260  YPTAU(I)  «  G 

DO  300  I  -  .1 1  NUP 

Q  *  0.0 

DO  301  J  ~  1,NC 

301  G  ~  Q  f  HC(.Tf J)*XC2( J) 

300  XWl(I)  -  G 

DO  302  I  t,NUP 

G  -  0.0 

DO  303  J  ~  1 » NOP 

303  G  «  G  +  ECdf J)*YPTAU<J) 

302  XW2  < I )  »  G 

DO  304  I  a  1 t NUP 

304  YC2(  I )  =  XWKI)  +  XU2<I) 

DO  305  I  a  1,NC 

Q  »  0.0 

DO  306  J  a  1  f  NC 
306  Q  =  G  +  FC(I> J)*XC2<J) 

305  XWKI)  ®  q 

DO  307  I  a  1 » NC 
G  a  0.0 

DO  308  J  =  :l»  NOP 


308 

G  =  G 

+ 

GC< I f J)#YPTAU< J) 

307 

XW2( I ) 

G 

DO  309 

I 

-  1  f  NC 

309 

XC2  < I ) 

r= 

XWKI) 

+  XW2 ( I ) 

DO  290 

I 

«  If  NUP 

290 

EKI)  i 

3 

YC3(  I )  - 

YC2 ( I ) 

DO  432 

I 

-  If NUP 

432 

YP2 ( 1 ) 

i= 

YP(  I )  - 

YPKI) 

DO  405 

I 

*  IfNP 

Q  -  0.0 

DO  406  J  ~  1 > NP 

4 06  Q  «  G  +  PHI  T  <  I  r  J )  #XP1  <  J  ) 

405  XWKI)  =»  Q 

DO  407  I  =  l.NP 

Q  «  0.0 

DO  408  J  »  If NUP 
408  G  a  Q  +  PSIT(I»J) #UP ( J ) 

407  XW2< I >  a  Q 

DO  409  I  =  KNP 

a  »  o.o 
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HO  410  J  =  lfNUP 

410  Q  =  Q  +  PSIT2<If J)*YP2< J> 

409  XW3< I )  =  Q 

HO  411  I  =  IfNP 

411  XPl(I)  =  XW1 < I >  +  XU2( 1 )  +  XW3 ( I ) 
no  433  I  =  if  NOP 

Q  »  0,0 

no  434  J  -  IfNP 
434  Q  =  Q  +  CP< I  * J)*XP1 < J) 

433  YPKI)  *  Q 

373  no  500  I  **  IfNP 

Q  «  0.0 

no  501  J  =  IfNP 

501  Q  =  Q  +  PHI T <  I»  J)#XP<  J) 

500  XWl(I)  -  0 

no  502  I  ~  IfNP 

Q  =  0.0 

no  503  J  =  lfNUP 

503  Q  =  0  +  PSIT  < I f J)#UP<  J) 

502  XU2< I )  *  0 

no  504  I  ~  IfNP 

504  XP(I)  =  XWl(I)  +  XU2< I ) 
no  507  I  --  If  NOP 

Q  =  0.0 

no  508  J  =  IfNP 
508  Q  =  Q  +  CP< If J)*XP( J) 

507  YP(I)  =  Q 

no  700  I  *  lfNUP 

Q  -  0.0 

no  701  J  *  1 ,NC 

701  Q  =  Q  +  HC<  I  f  J )  KtXCl  <  J ) 

700  XWl(I)  »  Q 

no  702  I  ■  lfNUP 

Q  -  0.0 

DO  703  J  »  If NOP 

703  Q  ■  Q  +  EC< If J)*YP< J) 

702  XW2 ( I )  *  Q 

DO  704  I  *  lfNUP 

704  YCl(I)  =  XWl(I)  +  XU2( I ) 

DO  705  I  -  IfNC 

Q  =  0.0 

DO  706  J  *  IfNC 

706  Q  «  Q  +  FC( I f J)*XC1 < J) 

705  XW1<I>  -  Q 

DO  707  I  =  IfNC 

Q  *  0.0 

DO  708  J  ■  If NOP 

708  Q  -  Q  +  GC<If J>*YP(J> 

707  XU2(  X  >  Q 

DO  709  I  =  IfNC 

709  XC1<I)  *  XWl(I)  +  XW2( X ) 

NN1  ■  NN1  +  1 

DO  221  I  *  If NOP 
221  UP<I>  ■  W<NN1 )~YC1 < I ) 
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0  *  0.0 

DO  521  J  ~  1 »NP 

521  0  =  0  +  PHTAUKIf  J)*XP1<  J) 

520  XW1 ( I )  =  Q 

DO  522  I  =  1 »NP 

0  *  0,0 

DO  523  J  *  1 t NUP 

523  0  =  0  +  PSTAUKIf  J)*UF‘<J) 

522  XW2 ( I )  =  0 

DO  524  I  =  lrNP 

524  XPTAU1 ( I )  =  XW1 < I )  +  XW2(I) 
DO  525  1  =  1 »NOP 

0  =  0.0 

DO  526  J  =  1 t NP 

526  0  =  0  +  CP< I f J)#XPTAU1 <  J) 

525  YPTAU1 <  I )  ■  0 

DO  529  I  =  1 t NUP 

0  =  0.0 

DO  530  J  -  1 »NC 

530  0=0+  HC(I, J)*XC3< J) 

529  XWKI)  =  Q 

DO  531  I  =  1 t NUP 
0  =  0.0 

DO  532  J  =  1 f NOP 

532  0=0+  EC<If J)*YPTAU1(J) 

531  XU2( I )  =  0 

DO  533  I  =  1 f  NUP 

533  YC3<  I  >  =  XWKI)  +  XW2(I) 

DO  535  I  =  1 rNC 

0  a  0.0 

DO  536  J  =  1 f NC 

536  0=0+  FC(If J)*XC3< J) 

535  XWKI)  »  0 

DO  537  I  =  1  f  NC 

Q  =  0,0 

DO  538  J  «  1 f NOP 

538  0  =  0  +  GCdf  J)*YPTAU1(  J) 

537  XW2( I )  «  0 

DO  539  I  =  1,NC 

539  XC3< I )  =  XW1<I)  +  XW2 ( I ) 

DO  540  I  =  lr NUP 

540  E2< I )  =  YC3< I )  -  YC2<I) 

DO  560  I  =  If NUP 

560  EA<  III )  =  EKI) 

220  CONTINUE 

SMEANA  =0,0 
DO  561  I  =  51 f 60 

561  SHE ANA  ■  SHEANA  +  EA<I) 

SHEANA  »  SMEANA/ 10 

VEA  =  0.0 
DO  563  I  =  51 f  60 

563  VEA  *  VEA  +  <EA< I ) -SMEANA ) #*2 
P£ASS( XCOUNT )  =  VEA/ 10 
RETURN 
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END 

SUBROUTINE  RANDUC IX » IY» YFL ) 
IY  =  IX#65539 
IF<IY>5»6,6 

IY  «  IY  +  2147483647+1 
YFL  »  IY 

YFL  =  YFL*0 . 4656613E-9 

RETURN 

END 
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EXAMPLE - - -2rH  ORDER  PLANT t  1ST  ORDER  CONTROLLER 

NO.  OF  PLANT  STATES  =  2 

NO.  OF  PLANT  INPUTS  «  1 
NO.  OF  EXTERNAL  INPUT  =■=  1 
NO.  OF  PLANT  OUTPUTS  =  1 

NO.  OF  CONTROLLER  STATES  (  EACH  CONTROLLER '  =  1 

PLANT  STATE  MATRIX  --  AP 

.0  1.0 

.0  -10.0 

PLANT  CONTROL  INPUT  MATRIX  —  BP 

.0 

200.0 


PLANT  OUTPUT  MATRIX  CP 

1.0  .0 


CONTROLLER  STATE  MATRIX  --  FC 
.523810 

OBSERVER  MATRIX  --  GE 
20.0 
25.0 

CONTROLLER  CONTROL  INPUT  MATRIX  •—  GC 
-.18162 

CONTROLLER  OUTPUT  MATRIX  (STATES)  HC 

1.0 

CONTROLLER  OUTPUT  MATRIX  < INPUTS)  —  EC 
1.381 

NT  -  51 

T  ■-  SAMPLE  PERIOD  -  0.0125  SEC 

DELTA  »  T / ( NT- 1 )  »  INCREMENT  USED  IN  .THE  NUMERICAL 

INTEGRATIONS  TO  COMPUTE  PSITAU.PSI  i',PS  I72r 
PSITAIJI  USING  TRAPEZOIDAL  RULE, 

W  IS  THE  EXTERNAL  INPUT  (WHITE  GAUSSIAN  NOISE  WITH 
MEAN  =  0,0.  AND  VARIANCE  «  1.0 

THE  SKEW  OF  THE  SECOND  CONTROLLER  IS  0,00725 

THE  NUMBER  OF  ITERATIONS  IS  5 


m 
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